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ABSTRACT 

New  concepts  in  aerospace  travel  have  renewed  interest 
in  modeling  and  compensating  for  the  effects  of  upper 
atmospheric  drag.  In  particular,  the  SDI  constellation 
requires  strict  orbital  element  maintenance.  This  thesis  is 
a  qualitative  verification  of  a  propellant  longevity  model 
for  low-altitude  earth  orbit  satellites  doing  intrack  micro- 
thrusting  to  overcome  atmospheric  drag.  The  original  model 
was  developed  at  the  Naval  Surface  Weapons  Center. 
Pertinent  orbital  mechanics  and  atmospheric  concepts  are 
reviewed.  The  model  and  its  computer  program  are  described. 
The  results  of  trend  and  sensitivity  analysis  reasonableness 
tests  are  presented.  Finally  suggestions  for  use  are  made. 
Many  plots  of  mission  life  predictions  are  presented.  The 
model  computer  program  and  sample  input  and  output  are  also 
included. 
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I.   INTRODUCTION 

All  current  low-earth-orbit  satellites  are  simply  placed 
in  orbit  and  then  tracked  as  their  orbits  decay  deeper  into 
the  atmosphere.  This  orbital  decay  is  due  to  forces  on  the 
satellite  produced  by  atmospheric  drag  and  the  non- 
spherical  earth. 

Modified  exponential-type  models  for  the  atmosphere 
produce  sufficiently  accurate  drag  predictions  to  track  the 
satellites. 

New  concepts  in  scientific  and  military  uses  of  space 
require  greater  accuracy  in  predicting  atmospheric  drag 
effects  then  produced  by  the  modified  exponential  models. 
These  new  uses  include  the  Aero-assisted  Orbital  Transfer 
Vehicle  (AOTV) ,  orbiting  space  stations,  the  Trans 
Atmospheric  Vehicle  (TAV) ,  and  the  strategic  defense 
initiative  satellite  constellation. 

The  environment  of  all  of  these  vehicles  is  called  the 
thermosphere ,  an  atmospheric  layer  studied  in  depth  only 
recently.  An  accurate  empirical  model  of  the  thermosphere 
has  also  been  produced. 

One  idea  under  investigation  at  the  Naval  Surface 
Weapons  Center  is  the  placement  of  many  low-earth-orbital 
defensive  satellites  in  precisely  controlled  and  maintained 
orbits.    This  set  of  defensive  satellites  is  termed  a 


constellation.  All  low  earth-orbital-satellites  drift  from 
their  initial  orbits  due  to  perturbating  accelerations 
caused  primarily  by  the  aspherical  earth  and  atmospheric 
drag.  In  order  to  be  functional,  the  orbits  of  the 
satellites  in  an  SDI  constellation  must  be  exactly 
maintained.  The  drift  of  a  satellite  must  be  cancelled  out 
by  using  an  onboard  propulsion  system  to  produce  forces  to 
offset  the  perturbating  accelerations. 

One  model,  proposed  by  Dr.  Dan  Parks  of  the  Naval 
Surface  Weapons  Center,  calculates  the  mission  lifetime 
expected  of  a  satellite  doing  intrack  micro-thrusting  to 
overcome  atmospheric  drag  effects.  The  end  of  mission  life 
corresponds  to  depletion  of  all  the  maneuvering  fuel.  This 
model  is  here-after  referred  to  as  the  propellant  longevity 
model.  The  propellant  longevity  model  incorporates  a  Bessel 
function  expansion  for  the  fuel-mass  decrement  on  each 
revolution.  However,  the  model  does  not  calculate  the  fuel 
necessary  to  cancel  out  orbital  drift  due  to  the  aspherical 
earth. 

This  thesis  gives  a  verification  of  the  low-altitude 
earth  satellite  propellant  longevity  prediction  model 
developed  by  Dr.  Parks.  The  theory  for  the  model 
development  itself  is  presented  in  Naval  Surface  Weapons 
Center  Technical  Report  (NSWC  TR)  83-243  (Parks  [Ref.  1]). 
Computer  code  for  the  model  was  developed  at  NSWC  under  the 
direction  of  Dr.  Parks.   In  addition  to  the  computer  code  of 


the  model,  several  other  NSWC  publications  with  applications 
to  the  model  were  forwarded  to  the  Naval  Postgraduate 
School. 

The  Thesis  is  organized  as  follows.  Chapter  I  presents 
a  general  introduction  to  the  propellant  longevity  model  and 
discusses  the  changes  made  to  the  computer  code  for  use  at 
the  Naval  Postgraduate  School.  A  basic  description  of  the 
orbital  mechanics  necessary  to  operate  the  model  is  provided 
in  Chapter  II.  Chapter  II  also  presents  an  overview  of  the 
atmospheric  physics  necessary  to  understand  the  density  and 
compositional  variations  in  the  thermospheric  layer  of  the 
thermosphere.  A  description  of  the  propellant  longevity 
model  is  given  in  Chapter  III.  Instructions  for  operation  of 
the  computer  code  are  also  included  there.  Chapter  IV 
contains  the  results  of  the  verification  process  sensitivity 
analysis  and  reasonableness  tests.  A  summary  is  given  in 
Chapter  V.  Appendix  A  contains  the  propellant  longevity 
Fortran  computer  code  for  the  IBM  3033.  Appendix  B  shows 
input  sets  and  their  corresponding  outputs.  Graphical 
presentation  of  the  tabular  results  is  included. 

The  coding  of  the  model  by  NSWC  was  completed  in  Fortran 
for  the  Control  Data  Corporation  (CDC)  6700  computer. 
Because  of  software  and  hardware  differences,  extensive 
modifications  to  the  model  coding  were  required  to  install 
the  program  as  an  operational  model  on  the  IBM  303  3  computer 
at  the  Naval  Postgraduate  School.   To  optimize  the  utility 


of  the  model,  some  adjustments  were  also  made  to  the 
original  computer  coding  of  the  propellant  longevity  model. 
Great  care  was  taken  when  modifying  the  code  to  preserve  the 
model  theory  as  contained  in  NSWC  TR  83-243. 

The  alterations  to  the  computer  program  fall  into  two 
categories:  changes  to  transport  and  implement  the  model  on 
our  IBM  3033,  and  changes  made  to  mainstream  the.  model 
verification  process.  Changes  to  transport  the  model  from 
the  CDC  environment  to  the  IBM  environment  involved  some 
logical  revision  and  syntax,  hardware  word-size 
considerations . 

The  hardware  word-size  problems  were  extensive.  The  CDC 
6700  uses  a  60-bit  word  size  for  single  precision 
calculations,  whereas  IBM's  standard  is  a  30-bit  word.  The 
propellant  longevity  model  requires  60  bit  precision.  Thus 
all  real  numbers,  library  functions,  and  input  constants 
were  promoted  to  double  precision  accuracy.  Another 
precision-related  problem  involved  the  use  of  fractional 
exponentiation.  All  one-half  exponents  were  changed  to  call 
the  SQRT  function  to  increase  precision  and  speed. 

Syntactical  changes  were  mostly  related  to  the  software 
environment  transportation  from  CDC  Fortran  to  IBM  Fortran 
77.  Library  functions  calls,  and  input  and  output  format 
statement  changes,  were  the  bulk  of  the  syntactical  changes 
made.  (Examples  of  library  function  changes  are  the  ARQTAN 
function  changed  to  the  ATAN2  function,  the  BESI  Bessel 
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function  call  changed  to  the  MMBSIR  Bessel  function  call  and 
so  forth) . 

Examples  of  changes  in  the  logic  to  in  order  to 
transport  and  implement  the  model  code  include  using  the 
Fortran  77  block  IF-THEN-ELSE  statement  for  the  CDC  Fortran 
IF  statement  and  the  use  of  DIMENSION  statements  in  each 
subroutine  that  is  processing  a  particular  array. 

The  second  category  of  changes  to  the  program  addressed 
mainstreaming  the  model  verification  process.  The  model 
code  as  down-loaded  from  the  NSWC  tape  contained  little 
code-internal  documentation.  Looking  towards  future 
implementation  and  maintenance  of  the  model  code,  extensive 
documentation  was  added  as  code  internal  comments.  Most  of 
this  documentation  was  based  on  the  research  literature  of 
NSWC.  Extensive  reference  to  the  current  publications  is 
made  in  the  code  documentation  itself. 

Two  changes  were  made  in  the  way  the  code  solves  the 
model  in  order  to  implement  the  theory  with  greater 
precision.  The  first  change  turned  on  the  subroutine  THRST 
which  was  designed  to  allow  for  preplanned  changes  to  the 
semi-major  axis  of  the  orbit  at  the  expense  of  onboard 
maneuvering  fuel.  (The  subroutine  code  down-loaded  from 
tape  did  not  return  the  fuel  consumed  for  the  orbit  adjust, 
record  the  resultant  change  in  the  semi-major  axis,  or 
receive  the  correct  intended  adjust  parameter.  All  these 
functions  are  performed  in  the  IBM  3033  version.) 


A  second  change  in  the  code  alters  the  way  in  which  the 
satellite  velocity  relative  to  the  atmosphere  is  computed. 
In  the  original  code  this  parameter  had  to  be  entered  as 
part  of  the  input  file  and  was  assumed  to  be  constant. 
Since  this  relative  velocity  is  actually  a  function  of  other 
transient  orbital  parameters,  recomputation  is  now 
automatically  accomplished  for  each  satellite  revolution 
based  on  current  values  of  the  changing  orbital  parameters. 
In  place  of  the  satellite-to-atmospheric-relative-velocity 
input,  the  user  now  inputs  the  atmospheric-to-earth  relative 
rotation  rates.  The  physics  associated  with  this  change  is 
discussed  later  on,  in  the  section  describing  the 
atmospheric  model. 

A  few  other  changes  to  the  propellant  longevity  code 
customize  the  output,  or  state  the  input  element  initial 
conditions.  Several  control  inputs  now  allow  for  a  range  of 
model  predictions.  These  are  explained  in  Chapter  III. 

The  actual  model  verification  is  accomplished  in  three 
steps.  First  does  the  model  address  the  original  problem? 
Second,  does  the  model  meet  the  criteria  of  common  sense? 
Finally,  how  well  does  the  model  solution  compare  when 
tested  against  real-world  data? 

No  real-world  data  sets  exist  for  testing  the  propellant 
longevity  model.  The  only  low-earth-orbit  satellite  having 
a  propulsion  system  to  maintain  orbital  parameters  is  highly 
classified  under  the  auspices  of  the  United  States  Air 
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Force.  The  data  from  this  program  is  not  available  in  an 
unclassified  state.  In  lieu  of  actual  data,  several 
additions  to  the  program  code  were  made  to  assist  in  the 
model  verification  process.  The  major  adjustment  is  the 
addition  of  a  feature  to  allow  for  a  sensitivity  analysis. 
The  subroutine  SNALYS,  allows  the  user  to  linearly  scan  any 
combination  of  the  20  input  elements,  from  a  user  specified 
minimum  to  any  maximum,  with  user  specified  granularity. 
The  output  of  the  sensitivity  analysis  can  be  descriptive 
and  tabular,  or  graphical  in  nature.  The  SNALYS  subroutine 
allows  for  a  plot  of  model  prediction  against  the  entire 
range  of  the  input  elements.  Results  of  input  element 
sensitivity  analysis  are  presented  in  the  Chapter  IV.  The 
sensitivity  analysis  feature  allows  for  optimum  orbital- 
element  and  satellite-ballistic-coefficient  design. 

Thought  was  given  to  possible  model  installation  on 
other  computers.  The  model  is  large  and  requires  about  20 
minutes  of  dedicated  central  processor  time  to  complete  an 
average  run.  As  such  the  model  is  rather  limited  to 
mainframe  computers.  Initial  calculations  based  on 
processor  time  indicate  an  IBM  AT  microcomputer  with  the 
requisite  80287  co-processor  and  additional  RAM  memory  might 
require  about  ten  hours  to  complete  an  average  run. 

This  chapter  presented  an  overview  and  motivation  for 
the  propellant  longevity  .model.  Changes  made  to  the 
computer  program  at  the  Naval  Postgraduate  School  were  also 
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discussed.  The  next  chapter  presents  the  necessary 
background  in  orbital  mechanics  and  atmospheric  concepts  to 
operate  and  understand  the  model. 
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II.   BACKGROUND 

This  chapter  presents  the  orbital  mechanics  and 
atmospheric  concepts  required  for  the  understanding  and 
operation  of  the  propellant  longevity  model.  The  orbital 
mechanics  presented  in  the  first  section  discusses  only  the 
orbital  element  sets,  concepts,  and  coordinate  systems 
specific  to  the  model.  In  the  second  section  of  this 
chapter,  Thermosphere  specific  atmospheric  concepts  and 
cycles  are  presented. 

A.   ORBITAL  MECHANICS 

This  section  covers  the  orbital  mechanics  for  operating 
the  mass  decrement  model.  If  the  coordinate  and  orbital 
element  systems  discussed  are  familiar,  this  basic  orbital 
mechanics  section  can  be  skipped.  Orbital  mechanics 
references  are  cited  in  the  text. 

Early  celestial  mechanics  models  were  based  on 
increasingly  accurate  observations.  As  observations  and 
record  keeping  improved  so  did  the  empirical  model  fit  to 
the  observations.  The  model  that  fathered  modern  orbital 
mechanics  was  the  empirical  model  derived  by  Johannes 
Kepler,  1571-1630.  In  1619  Kepler  formulated  three  laws 
(models)  of  planetary  motion.  These  models  were  based  on 
life  long  observations  recorded  by  Tycho  Brahe,  1564-1601. 
Kepler's  models  were  formulated  as  follows: 

13 


1.  The  orbit  of  each  planet  is  an  ellipse  with  the  sun  at 
one  focus. 

2.  The  line  joining  a  planet  to  the  sun  sweeps  out  equal 
areas  in  equal  times. 

3 .  The  square  of  the  period  of  a  planet  is  proportional 
to  the  cube  of  its  mean  distance  from  the  sun. 

Later,  Sir  Isaac  Newton,  1643-1727  validated  Kepler's 
laws  by  establishing  their  equivalence  with  his  law  of 
universal  gravitational  attraction. 

Newton's  model  for  gravitational  attraction  states  that 
two  bodies  attract  with  a  force  that  is  proportional  to  the 
product  of  their  masses  and  inversely  proportional  to  the 
square  of  the  distance  between  them.  Symbolically,  the 
force  on  the  mass  m  due  to  the  mass  M  is  given  by: 


Fg   =   -  «|^  -    (1) 

r 


where  G  is  the  universal  gravitational  constant  and  r  is  the 
vector  from  the  center  of  mass  of  M  to  the  center  of  mass  m. 
r  is  the  magnitude  of  the  distance  vector  r. 

A  body  experiences  an  acceleration  due  to  the  resultant 
of  all  the  forces  acting  on  the  body.  If  only  a 
gravitational  force  acts  on  the  body,  then  the  acceleration 
of  the  body  is  given  by: 


G(J^+m)  - 
=   a   = —^ —  r  (2) 

r 
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In  this  model,  r"  is  the  acceleration  vector  on  m  caused  by 
gravitational  attraction  of  M  on  m.  The  remaining  variables 
are  defined  as  in  the  previous  model.  If  M  is  much  greater 
then  m  (as  in  an  earth  satellite  system)  ,  G(M+m)  is  very 
nearly  GM.  Another  constant  can  then  be  defined  called  the 
gravitational  potential  parameter  u: 

u  =   Gm  (3) 

The  gravitational  acceleration  of  a  satellite  in  earth 
orbit  can  then  be  expressed  as: 

r"  +  u/r^  *  r  =  0  (4) 

The  gravitational  field  is  conservative.  This  means 
that  an  energy  constant  of  motion  can  be  defined  for  a  body 
moving  in  a  gravitational  field.  For  a  two  body,  an- 
atmospheric,  elliptic  orbit  system  the  specific  mechanical 
energy  is: 

C  =  v^/2   -   u/r  (5) 

In  this  model  ^  is  the  specific  mechanical  energy,  v-^-  is 
the   satellite   tangential   speed,   and   u   and   r  are   as 

previously  defined.   Note  that  E  is  the  sum  of  the  kinetic 
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energy  per  unit  mass  of  the  satellite  and  its  potential 
energy  per  unit  mass. 

The  model,  (5) ,  requires  satellites  at  specific  radii  to 
have  specific  tangential  speeds.  As  r  increases,  v  must 
decrease  to  maintain  ^  constant.  This  orbital  energy  model 
is  important  when  modelling  the  effects  of  atmospheric  drag 
which  is  the  major  topic  of  the  next  section. 

Newton's  inverse  square  law  of  gravitational  force 
requires  gravitational  field  orbits  to  be  elliptical  in  the 
shape.  The  mathematics  of  the  ellipse  is  important  in  the 
description  of  orbits  and  Figure  1  presents  some  of  the 
important  parameters. 

The  parameters  shown  in  Figure  1  are  defined  as  follows: 

a:   the  semi-major  axis. 

b:   the  semi-minor  axis. 

e:   the  eccentricity  where 

e  =   SQRT  [1  -  (b/a)2] . 

C:   the  ellipse  center. 

O:   the  ellipse  focus  (earth  center) . 

p:   the  semi-latus  rectum. 

r:   the  distance  from  a  position  on  the  ellipse 
the  focus  0. 

v:   the  true  anomaly. 

rp:   the  radius  at  perigee  (in  an  earth  system) . 

ra;   the  radius  at  apogee  (in  an  earth  system) . 
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Figure  1.   The  Elliptical  Orbit 

A:   the  apogee  (point  farthest  from  the  focus) . 

P:   the  perigee  (point  closest  to  the  focus. 

No  real  orbits  occur  in  just  two  dimensions.  To  describe 
real  orbits  a  third  dimension  must  be  added.  A  suitable 
reference  must  also  be  chosen.  One  of  the  most  convenient 
reference  frames  is  the  geocentric  reference  depicted  in 
Figure  2 .   In  this  system  a  fixed  nonrotating  earth   forms 
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the  center  of  the  celestial  sphere.  The  earth's  projected 
north  pole,  south  pole,  and  equator  become  the  north 
celestial  pole  (NCP) ,  the  south  celestial  pole  (SCP) ,  and 
celestial  equator  respectively. 


SCP 


Figure  2.   The  Geocentric  Celestial  Sphere 

To  account  for  the  apparent  motion  of  the  celestial 
sphere  as  viewed  from  the  earth,  a  reference  point  on  the 
celestial  sphere  must  be  chosen.  Traditionally  this  has 
been  the  First  Point  of  Aries,  describing  the  constellation 
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where  the  vernal  equinox  occurred  when  the  celestial 
reference  was  defined.  The  vernal  equinox  does  not  now 
occur  in  the  constellation  Aries.  Tradition  still  uses  the 
first  point  of  Aries  as  the  Celestial  Meridian  (CM)  of 
reference.  Figure  3  illustrates  the  vernal  equinox 
reference.  The  vernal  equinox  is  defined  as  the  time  when 
the  sun  is  directly  over  the  equator  at  noon  in  the  spring. 


VERNAL 
EQUINOX 


WINTER 
SOLSTICE 


FALL 
EQUINOX 


Figure  3 .   The  Vernal  Equinox 
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Figure  4  is  a  representation  of  an  orbit  in  real  space 
using  the  geocentric  system.  The  points  where  the 
projection  of  the  satellites  orbit  on  the  celestial  sphere 
crosses  the  celestial  equator  are  called  the  nodes  of  the 
orbit.  The  ascending  node  (N)  is  where  the  projected  orbit 
crosses  the  celestial  equator  heading  for  the  north 
celestial  pole.  The  descending  node  is  where  the  projected 
orbit  crosses  the  celestial  equator  heading  for  the  south 
celestial  pole.  The  angular  measurement  from  the  reference 
celestial  meridian  to  the  ascending  node  along  the  celestial 
equator  is  called  the  longitude  of  the  ascending  node^  and 
usually  given  the  symbol  Q,  .  The  angular  measurement  from 
the  ascending  node  to  the  point  of  perigee  (P)  along  the 
projected  orbit  is  called  the  argument  of  perigee, 
symbolized  by  oj  .  Both  the  longitude  of  the  ascending  node 
and  the  argument  of  perigee  range  from  0  to  360  degrees. 
The  angle  from  the  plane  of  the  celestial  equator  to  the 
plane  of  the  projected  orbit  is  called  the  inclination  and 
denoted  by  i.  Values  for  the  inclination  range  between  0 
and  180  degrees. 

Inclinations  of  0  and  180  degrees  are  termed  equatorial. 
In  equatorial  orbits  the  longitude  of  the  ascending  node  is 
not  defined.  Orbits  with  an  inclination  of  90  degrees  are 
called  polar.  Prograde  "orbits  are  those  with  an  inclination 
of  between  0  and  90  degrees  and  retrograde  orbits  have  an 
inclination  range  of  90  to  180  degrees. 
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NCP 


Figure  4.   A  Geocentric  Orbital  System 

Historically,  it  was  difficult  to  tell  much  about  the 
orbit  specific  parameters  from  the  earth  vantage  point. 
However,  a  value  that  could  be  measured  was  the  time  for  one 
complete  orbit.  If  the  radius  vector  moved  uniformly  along 
the  projected  orbit  with  the  average  rate  of  mean  motion  n 
(n  is  also  referred  to  as  the  mean  angular  velocity)  at  an 
epoch  T,   it  would  be  at  a  position  M  called  the  mean 
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anomaly.  The  mean  anomaly  has  a  range  of  0  to  3  60  degrees. 
It  is  measured  from  the  perigee  along  the  projected  orbit. 
The  mean  anomaly  corresponds  to  positions  on  the  orbital 
ellipse.  An  eccentric  anomaly  (E) ,  is  defined  as  the  angle 
between  the  projected  position  and  the  point  of  perigee 
subtending  at  the  center  of  the  ellipse.  Figure  5 
illustrates  eccentric  anomaly. 


Figure  5.   Orbital  Anomalies 
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The  difference  between  the  eccentric  anomaly  and  the  mean 
anomaly  is  important.  The  mean  anomaly  is  a  mathematical 
concept  corresponding  to  an  imaginary  position  of  a 
satellite  past  the  perigee  if  it  moved  at  an  average  angular 
velocity.  The  eccentric  anomaly  is  the  actual  projected 
angular  position  of  a  real  satellite  on  the  celestial 
sphere.  The  eccentric  and  mean  anomalies  are  equal  in 
circular  orbits  (e  =  0)  .  The  eccentric  and  mean  anomalies 
are  related  by  Kepler's  equation: 

M  =   E  -  e*sin  E  (6) 

The  solution  to  model  (6)  has  been  the  source  of  much 
investigation.  By  a  solution  it  is '  meant  for  a  given  M, 
determine  e  and  E.  In  1817  Friedrich  Wilhelm  Bessel,  1784- 
1846  invented  Bessel  functions  to  provide  an  infinite  series 
approximation  for  the  solution  of  the  model  (6)  .  Newton 
developed  a  converging  algorithm  for  its  solution.  The  mass 
decrement  submodel  of  the  propellant  longevity  model  uses 
Newton's  method  for  the  solution  of  Kepler's  equation. 

All  anomalies  are  concerned  with  the  measurement  of  the 
position  of  a  body  in  orbit.  All  anomalies  are  measured 
angularly  from  the  point  of  perigee  in  the  orbital  plane, 
and  vanish  at  perigee. 

When  e  and  E  have  been  found  from  Kepler's  equation  (6), 
then  using  the  equations  for  an  ellipse  in  Taff  [Ref.  2:p. 
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67],  the  remaining  elliptical  parameters  can  be  solved  for 
when  any  one  of  a,  b,  rp,  or  ra  is  known. 

The  minimum  number  of  parameters  needed  to  describe  an 
orbit  is  called  an  orbital  element  set.  Five  parameters  are 
needed  to  describe  an  orbit.  One  more  is  needed  to  position 
an  object  in  the  orbit.  Many  orbital  element  sets  have  been 
used  and  are  in  use.   The  classical  orbital  elements  are: 

semi-major  axis  (a) 

eccentricity  (e) 

inclination  (i) 

longitude  of  the  ascending  node  (Q,) 

argument  of  perigee  (co) 

time  of  perigee  passage  or  epoch  (T) . 
The  Keplerian  set  replaces  the  time  of  perigee  passage  with 
the  mean  anomaly  (M)  .  The  orbital  element  set  used  in  the 
mass  decrement  submodel  is  the  Kozai  mean  element  set,  very 
similar  to  the  Keplerian  set.  The  name  "Kozai  mean  element 
set"  refers  to  the  Kozai  perturbation  approximations  used  to 
solve  for  the  Keplerian  orbital  element  set.  The 
perturbation  theory  approximation  of  Y.  Kozai  is  discussed 
later. 

Up  to  now  all  orbital  considerations  have  been  limited 
to  the  classical  two-body  problem.  The  two  body  problem 
limits  the  forces  on  the  orbiting  body  to  the  gravitational 
force.  The  two-body  problem  also  requires  all  masses  be 
concentrated  at  a  single  point.    A  real  earth  satellite 
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experiences  many  forces  resulting  in  many  discrete 
accelerations.  The  theory  of  perturbations  was  developed  to 
account  for  the  motion  of  real  satellites. 

In  perturbation  theory  the  original  orbit,  with  only  the 
point-mass  gravitational  acceleration,  is  called  the 
osculating  orbit.  Orbital  elements  describing  the 
osculating  orbit  are  the  osculating  element  set.  The 
osculating  orbit  is  then  adjusted  periodically  to  account 
for  the  action  of  all  the  nonpoint-mass  gravitational 
accelerations. 

Accelerations  result  from  forces.  Examples  of 
perturbating  forces  on  an  earth  satellite  are:  oblate  earth 
gravitational,  atmospheric  drag,  sun  gravitational,  moon 
gravitational,  solar  radiation  pressure,  and  satellite  motor 
thrust. 

To  remain  in  a  closed  orbit,  the  point-mass 
gravitational  force  must  be  the  dominant  force  on  the 
satellite.  For  low  earth  orbit  satellites,  the  oblate  earth 
gravitational  force  is  the  next  greatest  force  and  is  termed 
the  dominant  perturbation.  Satellites  below  600  kilometers 
in  altitude  experience  an  atmospheric  drag  force  of  the  same 
order  of  magnitude  as  the  oblate  earth  force.  The 
atmospheric  drag  and  oblate  earth  effects  are  so  large  for 
low  earth  orbit  satellites  that  the  effects  of  the  other 
perturbations  (except  motor  thrust)  are  .  lost  in  the 
background  error. 
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Central  to  the  problem  of  perturbation  theory  is  the 
that  the  three  body  problem  has  never  been  solved 
mathematically  in  general  closed  form.  Most  modern 
approximations  involve  power  series  or  Legendre  polynomial 
expansions.  The  state  of  the  art  in  orbit  perturbation 
prediction  only  uses  only  the  first-order  terms  of  these 
expansions,  termed  first-order  effects.  An  exception  is  a 
theory  advanced  in  1959  by  Y.  Kozai.  His  approximations  are 
computationally  efficient  and  allow  for  the  use  of  more  of 
the  oblate  earth  effects  (called  J2,  J3 ,  and  J4  terms  and 
discussed  later) .  Kozai ' s  theory  is  flawed  in  the  higher- 
order  terms  due  to  the  introduction  of  mathematical 
singularities.  Nevertheless,  proponents  for  its  use  argue 
that  it  does  work  and  predicts  the  oblate  earth  effects  to  a 
reasonable  degree  of  accuracy.  The  mass  decrement  submodel 
of  the  propellant  longevity  model  presented  in  this  thesis 
uses  Kozai 's  theory. 

The  gravitational  potential  parameter  u  was  presented  as 
singularly  determining  the  gravitational  potential  in  the 
earlier  development  of  the  gravitational  acceleration  model. 
In  earth  systems,  the  gravitational  potential  is  a  complex 
function  of  the  earth's  mass  distribution. 

The  earth  is  not  a  point  mass.  Neither  is  it  a 
homogenous  sphere.  It's  true  mass  distribution  is  complex. 
The  real  distribution  of  matter  on  the  earth  forms  a  body 
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termed  the  geoid.     This  geoid  determines  exactly  the 
gravitational  potential  function. 

The  determination  of  the  geoid  is  based  on  ground  and 
satellite  measurements.  Many  military  programs  depend  on 
the  precise  position  of  a  point  on  the  surface  of  the  earth. 
The  military  has  led  the  way  in  geoid  determination  in  the 
past  two  decades. 

Geoidal  studies  reveal  important  facts  about  the 
aspherical  shape  of  the  earth.  A  fluid  spinning  sphere 
should  exhibit  a  bulge  at  the  equator  and  flattening  at  the 
poles.  This  is  termed  oblateness.  The  oblate  deviation  is 
the  dominant  deviation  for  the  earth  as  a  sphere. 
Measurements  of  the  earth's  oblateness  indicate  that  it  is 
more  oblate  than  the  current  spin  rate  could  produce  by 
itself.  This  excess  oblateness  indicates  two  factors:  the 
earth's  core  is  more  solid  than  previously  thought,  and  the 
earth's  spin  has  slowed. 

The  earth  deviates  in  other  ways  from  being  a  sphere. 
Since  all  deviations  can  be  expressed  in  terms  of  sines  and 
cosines  they  are  referred  to  as  spherical  harmonics. 
Harmonics  are  divided  into  categories.  Zonal  harmonics  are 
symmetric  about  the  earth's  spin  axis  and  are  not  longitude 
dependent.  Sectorial  harmonics  are  only  longitude 
dependent.  Harmonics  dependent  on  both  longitude  and 
latitude  are  termed  tesseral.  Odd  numbered  harmonics  are 
antisymmetric  about  the  equatorial  plane;    even  numbered 
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harmonics  are  symmetric  about  the  equatorial  plane. 
Harmonics  are  also  classified  by  number.  The  number  refers 
to  the  subscript  of  the  J  coefficient  (defined  later) . 

Bate  et  al.,  [Ref.  3],  offers  a  simplified  development 
of  the  gravitational  potential  function.  If  the 
gravitational  potential  function  0  of  a  geoid  is  known,  the 
acceleration  due  to  this  gravitational  function  is  modeled 
by: 


:n   — 


=  V$  (7) 

One  common  expression  for  the  gravitational  potential  or 
geoid  function  $  is: 

OO  ]^ 

$  =  li[l  -  y  J  (_^)n  p  3i^  L]  (8) 

n=2  ^  ^    ^ 

In  this  model  u  is  the  gravitational  potential  parameter,  Jn 
is  an  experimentally  determined  J  coefficient,  r^  is  the 
radius  of  the  earth,  Pn  is  a  Legendre  polynomial  of  order  n, 
L  is  the  geocentric  latitude,  and  r  is  the  radius  of  the 
satellites  position.  Taking  the  gradient  of  the  potential 
function  yields  an  expression  for  the  acceleration  due  to 
the  gravitational  potential.  Bate  et  al.,  [Ref.  3:pp.  419- 
423]  details  the  resultant  potential  induced  acceleration  to 
the  first  seven  terms. 
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Recent  studies  have  determined  the  geoid  function  out 
through  the  J44  term.  Values  used  by  the  mass  decrement 
model  by  the  subroutine  GEOP  are  1082.76*10"^,  -2.55*10"^, 
and  1.56*10"^,  for  the  J2 ,  J3 ,  and  J4  terms  respectively. 
The  J2  and  J4  terms  are  even  numbered,  and  refer  to  even 
harmonics  symmetric  about  the  equator.  The  J3  term  is  odd 
numbered,  and  refers  to  odd  harmonics  which  are 
antisymmetric  about  the  equator.  The  J2,  J3 ,  and  J4  terms 
are  all  zonal. 

Current  perturbation  theory  expresses  $  to  the  J2 
term.  The  theory  advanced  by  Kozai  in  1959  appears  to  use 
the  J2,  J3,  and  J4  coefficients  in  a  first-order 
perturbation  theory  approximation  (because  singularities 
occur  in  terms  higher-order  than  J2) .  Taff  [Ref.  2:pp.  332- 
342]  takes  strong  exception  to  the  theory. 

The  subroutine  GEOP  in  the  propellant  longevity  model 
uses  Kozai 's  theory  to  calculate  the  geoidal  effects  on 
orbital  parameters  through  the  first  three  zonal 
coefficients  (J2,  J3,  and  J4  terms).  Their  effect  is 
calculated  once  each  orbital  revolution,  and  the  Kozai  mean 
element  set  is  adjusted  accordingly.  While  using  the  model, 
it  is  important  to  understand  that  as  the  altitude  of  the 
satellite  decreases,  The  aspherical  earth  perturbations 
become  much  greater.  The  orbital  element  drift  due  to  the 
aspherical  earth  is  thus  much  faster  in  lower  altitude 
orbits. 
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In  summary,  the  computer  program  of  the  propel lant 
longevity  model  uses  three  orbital  element  sets  in  a 
geocentric  reference.  The  input  orbital  element  set  is  the 
Brouwer  mean  element  set  of  the  Naval  Space  Surveillance 
System  (NAVSPASUR) ^  system.  Chapter  III  and  Appendix  A 
explain  the  input  elements  in  detail  in  the  input 
description.  The  Brouwer  mean  elements  are  first  converted 
into  a  Brouwer  osculating  set.  These  osculating  orbital 
elements  are  then  transformed  into  a  Keplerian-like  Kozai 
mean  set.  Half  of  the  computer  program  of  the  propellant 
longevity  model  transforms  the  Brouwer  mean  element  set  into 
the  Kozai  mean  set.  These  element  set  transformations  are 
used  to  prepare  real-world  orbital  element  sets  for  use  in 
the  propellant  longevity  predictions.  Once  the  computer 
program  has  solved  for  the  Kozai  mean  element  set  all 
perturbation  effects  result  in  adjustments  to  the  Kozai  mean 
set.   Adjustments  are  made  once  per  revolution. 

B.   THE  ATMOSPHERE  MODEL 

Central  to  the  validation  of  any  model  is  an 
understanding  of  the  system  to  be  modeled.  The  model  under 
investigation  determines  the  propellant  longevity  for  low 
earth  orbital  satellites  doing  in  track  micro-thrusting  to 


^The  NAVSPASUR  system  publishes  a  semiannual  report  on 
the  1000+  orbiting  objects  they  track  in  The  Satellite 
Situation  Summary.  The  orbits  are  described  in  a  five  card, 
80  column,  Fortran  style  input  format.  The  report  is 
available  on  request  from  the  Commanding  Officer,  Naval 
Space  Suirveillance  System,  Dahlgren,  Virginia. 
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overcome  the  effects  of  atmospheric  drag.  The  fuel  required 
by  a  particular  rocket  motor  on  a  satellite  to  perform  a 
Icnown  amount  of  work  can  be  computed  from  the  rocket 
equation.  The  problem  occurs  in  computing  the  amount  of 
work  required. 

To  compute  the  required  work  by  the  satellite  propulsion 
system  to  overcome  atmospheric  drag  requires  knowledge  of 
the  atmospheric  environment  of  the  satellite.  For  low  earth 
satellites  the  atmospheric  environment  is  the  thermosphere . 
The  physics  of  the  thermosphere  is  discussed  in  this 
section.  Before  discussing  the  thermosphere,  some 
preliminary  thermospheric  properties  are  presented. 

For  low  earth  orbital  satellites,  the  drag-induced 
change  in  the  orbital  elements  is  significant.  The  oblate 
earth  (referred  to  as  the  "J2  term")  orbital  perturbation  is 
normally  the  dominant  perturbation  term.  The  perturbations 
due  to  drag  are  of  the  same  order  of  magnitude  as  those  due 
to  the  oblate  earth  for  orbits  below  600  kilometers. 

The  drag  force  caused  by  the  satellite  moving  through 
the  atmosphere  always  acts  against  the  motion  of  the 
satellite.  Neglecting  the  lift  effects  due  to  the 
satellite's  shape  and  attitude,  the  drag  force  causes  work 
to  be  done  on  the  atmosphere  by  the  satellite. 

From  Bate  e  al.,  [Ref.  3:pp.  15-16]  a  satellite  in  orbit 
has  a  constant  specific  mechanical  energy  C  given   by 
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2 

r   =   It_  -  ii  (9) 

^       2     r 


where  B,  is  the  mechanical  energy  per  unit  mass,  v^  the 
satellite  speed,  u  the  gravitational  potential  constant,  and 
r  the  radius  of  the  orbit. 

The  above  model  says  that  the  mechanical  energy  per  unit 
mass  of  a  body  in  an  an-atmospheric,  two  body,  circular 
orbit  is  the  sum  of  its  kinetic  energy  per  unit  mass  and 
gravitational  potential  energy  per  unit  mass. 

Conservation  of  energy  requires  that  the  work  done  on 
the  atmosphere  by  the  satellite  must  reduce  the  orbital 
energy  of  the  satellite.  The  differential  element  of  work 
done  on  the  atmosphere  by  the  satellite  results  in  a 
differential  reduction  in  the  specific  kinetic  energy  of  the 
satellite.  Assuming  the  specific  mechanical  energy  is  a 
constant,  the  radius,  r,  of  the  orbit  must  therefore 
decrease.  As  the  radius  decreases,  the  velocity  tends  to 
increase  in  order  to  maintain  the  an-atmospheric 
conservation  of  specific  mechanical  energy. 

In  non-circular  (real  elliptical)  orbits,  the 
atmospheric  effects  are  greatest  at  perigee  (Figure  1)  . 
Orbital  mechanics  requires  a  differential  reduction  in 
orbital  velocity  at  perigee  to  result  in  a  reduction  in  the 
semi-major  axis  and  orbital  eccentricity.  The  major  effect 
of  transatmospheric  satellite  motion  is  to  reduce  the  semi- 
major  axis  and  the  eccentricity  of  the  elliptical  orbit. 
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Taff  [Ref.  2:pp.  352-353]  develops  the  mathematical 
equations  for  the  drag-induced  changes  in  the  semi-major 
axis,  eccentricity,  argument  of  perigee  (Figure  4) ,  and  mean 
anomaly  (Figure  5)  .  The  changes  in  the  argument  of  perigee 
and  the  mean  anomaly  are  shown  to  be  periodic.  The  net 
effect  on  the  orbit  due  to  atmospheric  drag  is  to 
circularize  the  orbit  closer  to  the  earth,  with  little 
change  to  the  other  orbital  elements. 

An  accepted  model  for  the  magnitude  of  the  drag  force  as 
presented  in  Taff  [Ref.  2: p.  352]  is: 

A  pC       2 

Fd   =     o    V   /  (10) 

"       2m    rel 

In  this  equation,  Fq  is  the  drag  force  magnitude,  A  the 
effective  satellite  cross-sectional  area,  C^  the  drag 
coefficient,  m  the  satellite  mass,  p  the  atmospheric 
density,  and  v^^^  the  velocity  of  the  satellite  relative  to 
the  atmosphere. 

The  upper  atmosphere  is  dynamic.  The  changes  that  must 
be  modeled  to  achieve  a  valid  prediction  for  the  drag  force 
are  those  that  effect  the  drag  force  factors.  The  factors 
in  the  drag  force  model  that  are  dependent  on  the  atmosphere 
are  density,  a  drag  coefficient  (dependent  on  satellite 
physical  and  fluid-flow  characteristics) ,  and  the  relative 
velocity  of  the  satellite  in  the  atmosphere. 
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To  model  the  atmospheric  effects  on  density,  drag 
coefficient,  and  relative  velocity,  some  background  in  upper 
atmosphere  physics  is  needed.  The  descriptive  presentation 
here  is  carefully  developed  to  be  relative  to  the  above 
three  drag  force  dependent  factors.  A  more  comprehensive 
treatment  of  atmospheric  physics  can  be  found  in  a  work  of 
Fleagle  [Ref.  4].  An  excellent  technical  treatment  of 
thermospheric  dynamics  can  be  found  in  the  works  of  Jacchia 
[Refs.  5,6,7]. 

Jacchia  formulated  the  most  widely  used  thermospheric 
model  currently  employed  by  professional  orbital  analysts 
(known  as  the  Jacchia  J60  model  or  simply  J60) .  His  latest 
comprehensive  empirical  thermospheric  model  (Jacchia  J77  or 
J77)  is  many  orbital  analysts'  preferred  refinement  for 
modeling  the  atmosphere. 

The  Jacchia  J77  model  is  highly  tabular,  uses 
considerable  computer  time,  and  requires  an  extensive  input 
file.  However,  these  negative  aspects  are  balanced  by  the 
orders-of-magnitude  increase  in  accuracy  of  density 
predictions.  Jacchia  [Ref.  7:p.  2]  presents  a  table  of 
residuals  for  predicted  versus  observed  densities  for 
various  satellites,  as  well  as  time  intervals  for  the 
Jacchia  J77  atmosphere  model.  All  density  residuals,  are 
much  smaller  than  those  produced  by  exponential  models  in 
the  same  time  intervals.  When  extremely  accurate  density 
predictions  are  required  in  some  submodel,  the  use  of  the 
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Jacchia  J77  atmosphere  model  may  well  be  worth  the 
additional  computer  time  and  expense. 

A  simple  model  of  the  atmosphere  considers  it  as  a 
homogeneous  mass  of  constant  density.  A  refinement  is  made 
by  modeling  the  density  as  decreasing  exponentially  with 
altitude.  The  exponential  model  predicts  pressure  and 
density  will  fall  one  tenth  for  every  ten  kilometer  increase 
in  altitude.  Prior  to  1960  this  exponential  model  was  used 
extensively  for  upper  air.  calculations  involving  meteor 
trails.  With  the  advent  of  low  earth  orbital  satellites,  a 
greater  degree  of  accuracy  in  predicting  atmospheric  effects 
is  required.  The  predicted  satellite  decays  from  the 
exponential  model  are  very  different  from  the  observations. 
The  Jacchia  J60  atmosphere  model  was  formulated  to  improve 
decay  predictions.  It  was  based  on  measurements  form  high 
altitude  balloons,  sounding  rockets,  and  ground  stations. 

The  atmosphere  is  divided  roughly  into  concentric 
shells,  each  with  specific  properties.  The  properties  used 
in  this  subdivision  are  after  Fleagle.   [Ref.  4] 

From  the  early  days  of  atmospheric  studies,  temperature 
and  moisture  have  been  recorded.  The  division  system 
described  by  Fleagle  [Ref.  4]  is  based  primarily  on  a  height 
versus  temperature  profile. 

The  height  versus  temperature  profile  is  graphed  in 
Figure  6.  The  temperatures  and  altitudes  are  those  of  the 
1975  U.S.   standard  atmosphere.   This  temperature  behavior 
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Figure  6.   Atmospheric  Temperature  Profile 

approximates  general  conditions  (large  transient  variations 
have  been  observed) .  Figure  6  also  shows  the  boundaries 
between  concentric  shells,  called  pauses.  The  following 
description  of  the  atmospheric  layers  is  presented  in 
Fleagle  [Ref.  4:pp:  79-84]. 

The  region  of  the  atmosphere  closest  to  the  earth  is  the 
troposphere  and  contains  eighty  percent  of  the  atmospheric 
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mass.  This  region  is  characterized  by  near  linearly 
decreasing  temperature  with  altitude.  The  dynamics  are  due 
to  earth  surface  energy  transitions.  All  of  the  earth's 
weather  occurs  in  this  region.  The  height  of  the 
troposphere  is  variable  and  depends  on  both  latitude  and 
season.  Along  the  equator  the  maximum  altitude  is  16  to  18 
kilometers.  At  the  poles  the  troposphere  is  highly 
seasonal,  being  nearly  absent  during  the  winter  months  and 
extending  to  eight  to  ten  kilometers  in  the  summer.  Capping 
the  troposphere  is  a  region  of  constant  temperature  with 
altitude,  called  the  tropopause. 

Next  in  the  series  of  atmospheric  layers  is  the 
stratosphere .  The  stratosphere  is  characterized  by  linearly 
increasing  temperature  with  height.  The  physical  dynamics 
are  due  mainly  to  the  absorption  of  solar  ultra-violet 
radiation  by  the  ozone.  The  maximum  altitude  of  the 
stratosphere  is  about  50  kilometers  where  it  is  capped  by  an 
area  of  local  maximum  temperature  termed  the  stratopause. 

Above  the  stratopause  is  the  mesosphere.  This 
atmospheric  layer  exhibits  linearly  decreasing  temperature 
with  height.  The  energy  physics  of  the  mesosphere  are  not 
well  understood.  Bulk  mixing  of  the  gasses  still  dominates 
diffusion  as  thermodynamic  processes.  Any  molecular 
dissociation  of  atmospheric  molecules  produced  by  radiation 
absorption  is  rapidly  consumed  by  chemical  recombination  of 
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these  dissociated  species.  The  mesosphere  is  topped  by  the 
turbopause  at  an  altitude  of  about  90  kilometers. 

Above  the  turbopause  the  gross  physics  of  the  atmosphere 
changes  drastically.  This  region  is  the  environment  of  all 
current  low  earth  orbit  satellites.  Termed  the 
thermosphere .  it  extends  from  the  turbopause  to  an  altitude 
of  1000  to  2500  kilometers.  Here  bulk  molecular  mixing  is 
dominated  by  diffusion.  Lighter  gasses  diffuse  upwards  and 
heavier  ones  remain  at  the  lower  altitudes.  Disassociation 
of  molecular  species  becomes  increasingly  important,  this 
results  in  an  increase  of  individual  atomic  species. 
Standard  flow  equations  (Navier-Stokes)  are  extremely 
difficult  to  apply,  and  physical  descriptions  are  more 
discrete  (quantum  mechanical) .  Flight  by  flow-supplied  lift 
is  not  possible.  The  standard  methods  of  computing  drag 
coefficients  do  not  apply.  The  ideal  gas  approximation  is 
also  not  valid  (P*V  =  n*R*T  does  not  apply)  .  Above  the 
turbopause,  density  is  an  increasing  function  of 
temperature.  Only  about  one  millionth  of  the  atmosphere 
remains  above  the  tropopause. 

The  stratification  of  gaseous  components  in  the 
thermosphere  is  significant.  The  drag  coefficient  has  been 
shown  to  be  composition  dependent.  Thus,  as  the  gaseous 
composition  changes  the  drag  coefficient  must  be  adjusted  to 
reflect  the  observed  drag-induced  orbital  decays  (Jacchia 
[Ref.  7]).   Figure  7  shows  the  dominant  atmospheric   species 
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Figure  7,   Atmospheric  Composition 

by  altitude  with  the  accepted  drag  coefficients.  The  drag 
coefficients  are  predictions  obtained  analyzing  the  drag 
coefficient  needed  to  produce  the  observed  drag  effects  on 
orbital  decay  for  a  satellite  at  a  particular  altitude. 
Thus  the  drag  coefficient  is  a  proportionality  constant  in 
the  drag  force  model. 


39 


For  altitudes  below  110  kilometers,  molecular  nitrogen 
is  the  main  atmospheric  gas.  Above  this  region  the 
dissociation  of  the  molecular  oxygen  by  ultra-violet  light 
absorption  causes  atomic  oxygen  to  become  increasingly  the 
main  gaseous  constituent.  At  an  altitude  of  500  kilometers, 
helium  replaces  atomic  oxygen  as  the  main  atmospheric  gas. 
At  altitudes  above  1500  kilometers  atomic  hydrogen  becomes 
the  dominant  gas  (see  Figure  7) . 

Traditional  research  has  assigned  a  drag  coefficient  of 
2.20  to  altitudes  between  200  and  400  kilometers.  As  the 
altitude  of  the  satellite  increases  above  500  kilometers, 
corresponding  to  satellite  movement  through  helium,  a  drag 
coefficient  of  at  least  3.00  must  be  assigned  to  maintain 
agreement  with  observation  (see  Figure  7) . 

Jacchia  [Ref.  7]  lists  the  observed  variational  changes 
in  density  and  composition  in  the  thermosphere .  They  are 
the  two  component  solar  activity  related,  semiannual, 
diurnal,  seasonal  latitudinal,  and  geomagnetic. 

The  two  component  solar  thermospheric  variations  are 
consequences  of  changes  in  solar  ultra-violet  (EUV)  flux. 
The  EUV  flux  changes  periodically  due  to  two  solar  cycles. 
These  are  the  27  day  rotational  or  disk  cycle  and  the  11 
year  sunspot  or  active  area  cycle, 

EUV  flux  is  not  measurable  from  the  ground  due  to  almost 
complete  absorption  in  the  atmosphere.  Empirical  evidence 
has  suggested  a  strong  parallel  of  the  EUV  flux  and  the 
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ground  measurable  10.7  centimeter  flux  (F10.7)  which  has 
been  used  in  some  successful  thermospheric  models.  The 
Jacchia  1960  model  used  the  10.7  centimeter  flux  as  an  input 
to  model  thermospheric  density  resulting  in  an  improvement 
of  decay  predictions  over  the  exponential  model.  (Indeed, 
it  is  only  recently  that  there  has  been  a  requirement  for 
greater  accuracy  in  drag  predictions.) 

The  thermosphere  responds  to  increased  EUV  flux  by 
absorbing  more  energy,  causing  increased  heating.  This 
increased  heating  causes  an  increase  in  temperature  and 
density.  The  thermospheric  response  to  changes  in  the  EUV 
flux  is  nonlinear.  Increases  in  EUV  flux  due  to  disk 
(rotational)  effects  produce  more  thermospheric  heating  then 
a  like  increase  in  flux  due  to  active  area  effects. 
However,  Jacchia  noticed  in  1973  that  the  disk  component  of 
the  EUV  flux  can  be  linearly  related  to  an  average  10.7 
centimeter  flux.  He  recommends  using  a  71  day  average 
corresponding  to  three  solar  rotations  [Ref.  7].  Tables  for 
the  10.7  centimeter  flux  have  been  developed. ^  The  Jacchia 
1977  atmosphere  incorporates  the  average  10.7  centimeter 
flux  as  an  input,  used  as  an  indicator  of  the  thermospheric 
heating  capacity  of  the  sun's  radiation. 


^Orr  [Ref.  8: Appendix  B]  lists  values  compiled  and 
predicted  by  the  NASA  Marshall  Space  Flight  Center.  Values 
range  from  a  minimum  of  73.3  to  a  maximum  of  242.9.  The 
units  for  the  10.7  centimeter  flux  are  10"22 
watts/meter^/Hz . 
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The  semiannual  temperature  and  density  variation  has 
been  modeled  in  many  ways.  Jacchia  [Ref.  7: pp.  47-50]  fits 
a  modified  harmonic  empirical  model  to  the  biannual, 
sinusoidal  cycle  of  observed  density.  Maxima  occur  in  April 
and  October.  The  Jacchia  1971  and  1977  models  for  the 
semiannual  density  variation  are  based  on  observations  of  12 
years  of  satellite  data.  The  observed  density  variations 
are  thought  to  be  a  consequence  of  the  earth-sun 
orientation.  The  mechanism  for  the  observed  effect  is 
unknown,  but  may  be  related  to  the  maximum  equatorial 
heating  at  the  equinoxes.  Equinoxes  (sun  at  maximum 
altitude  over  the  earth's  equator)  occur  in  late  March  and 
September.  Figure  8  illustrates  the  semiannual  bulge  in 
temperature  and  density  observed  just  after  the  equinoxes. 
Due  to  the  earth's  oblation,  the  bulk  mass  of  the  atmosphere 
is  concentrated  at  the  equator.  The  most  direct  exposure  of 
the  equatorial  atmosphere  occurs  at  the  equinoxes.  The 
temperature  and  density  thus  peak  shortly  after  periods  of 
maximum  heating. 

Seasonal-latitudinal  thermospheric  variations  are  again 
earth-sun  orientation  related.  Variations  in  density  and 
composition  are  dependent  on  the  season  and  latitude.  The 
Jacchia  1977  model  fits  modified  harmonic  empirical  models 
to  observed  data.  The  seasonal  latitudinal  variations 
interact  strongly  with  the  semiannual  variation,  complicat- 
ing the  modeling.   The  winter  helium  bulge  is  the  best  known 


42 


SPRING    MAXIMUM 
^    EQUITORIAL     BULGE 


FALL    MAXIMUM 
EQUITORIAL    BULGE 


Figure  8.   The  Semiannual  Atmospheric  Density  Bulge 

effect  of  this  variational  factor.  Figure  9  shows  the 
general  shape  change  of  the  temperature  and  density  bulge 
with  season.  Note  that  during  the  June  solstice  the  bulge 
is  centered  at  about  45  degrees  north  latitude  whereas 
during  the  December  solstice  the  bulge  is  centered  at  about 
45  degrees  south  latitude. 


43 


SUMMER  BULGE 


WINTER  BULGE 


Figure  9.   The  Seasonal-Latitudinal  Atmospheric 
Density  Bulge 


The  diurnal  bulge  is  shown  in  Figure  10.  It  is  caused 
by  the  daytime  heating  of  the  side  of  the  atmosphere  facing 
the  sun.  Lagging  the  actual  (most  direct)  solar  exposure  by 
about  four  hours,  the  heated  side  shows  a  marked  increase  in 
temperature  and  density.  Many  current  models  track  and 
predict  the  position  of  the  diurnal  bulge. 
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Figure  10.   The  Diurnal  Atmospheric  Density  Bulge 

Geomagnetic  heating-  of  the  thermosphere  is  not  well 
understood.  At  times  of  high  geomagnetic  activity 
(corresponding  to  an  index  of  7)  large  heating  effects  are 
produced.  The  combined  effects  of  geomagnetic  and  EUV 
heating  can  produce  changes  in  density  of  three  orders  of 
magnitude  at  600  kilometers  (Jacchia  [Ref.  5:p.  82]). 
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In  summary,  the  thermosphere  is  a  very  dynamic  region  of 
the  atmosphere.  Many  of  the  observed  changes  require 
further  research  to  complete  analytical  models.  The  current 
modified  exponential  empirical  models  can  be  orders  of 
magnitude  in  error  in  their  density  predictions.  At  higher 
earth  orbits  orders  of  magnitude  error  in  density  is  of 
little  consequence.-^ 

The  model  chosen  for  density  and  compositional  effects 
in  the  atmosphere  depends  on  the  importance  of  the  short- 
term  cyclic  effects.  Computations  for  orbits  of  months  in 
duration  and  less  than  600  kilometers  in  altitude  will  be 
subject  to  large  errors  using  the  simple  empirical  models 
(Jacchia  1966) .  For  these  short  missions  the  more  accurate 
Jacchia  1977  model  seems  more  appropriate.  (Performance  of 
the  Jacchia  1977  empirical  model  has  been  remarkable  with 
density  errors  rarely  exceeding  two  percent.)  Missions 
lasting  for  decades  may  be  able  to  use  the  computationally 
more  efficient  Jacchia  J60  model  to  predict  very  long  term 
drag  effects. 

The  one  drag  force  dependent  factor  not  yet  discussed  is 
the  satellite-to-atmosphere  relative  velocity.  This 
relative  motion  is  influenced  by  three  gross  factors:   local 


■^Above  800  kilometers  the  drag  term  in  the  orbital 
decay  model  used  by  the  Satellite  Control  Facility  in 
Sunnyvale,  California  is  very  small.  Its  main  function  is  a 
catch-all  for  errors  in  other  perturbating  terms. 
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winds,   satellite   tangential   velocity,   and   atmospheric 
rotation. 

Thermospheric  wind  has  not  been  measured.  Mesospheric 
wind  has  been  measured  based  on  the  observation  on  meteor 
trails.  Fleagle  [Ref.  4: p.  82]  reports  meospheric  winds 
exceeding  0.150  kilometers  per  second.  This  may  sound 
large,  but  must  be  compared  with  the  satellite's  velocity 
(which  is  on  the  order  of  10  kilometers  per  second) . 

Satellite  tangential  velocity  is  a  major  factor  in  the 
satellite-to-atmosphere  relative  velocity.  When  combined 
with  the  rotating  atmosphere  it  accounts  for  the  bulk  of 
relative  motion. 

Most  earth  orbit  calculations  simplify  computation  of 
perturbation-induced  changes  in  the  orbital  elements  by 
assximing  a  fixed  earth.  When  considering  drag  a  non- 
rotating  earth  assumption  can  lead  to  a  severe  error.  The 
earth  is  not  fixed.  The  atmosphere  is  pulled  along  with  its 
rotation,  and  lags  the  earth  rotation  by  a  small  fraction. 
Figure  11  shows  the  rotating  atmosphere  effect.  For  low 
inclination  (prograde)  orbits  the  atmosphere  is  actually 
following  the  satellite.  In  high  inclination  (retrograde) 
orbits  the  atmosphere  is  moving  against  the  satellite.  When 
comparing  a  typical  16  revolution  per  day  satellite 
tangential  velocity  in  a  prograde  atmosphere  to  the  same 
satellite  in  a  retrograde  atmosphere,  a  15  percent  variation 
in  relative  velocity  can  occur.   This  is  compared  to  the 
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Figure  11.   Orbital  Inclinations  in  a 
Rotating  Atmosphere 


unknown   (but  probably  less  than  two  percent  effect)   of 
atmospheric  winds. 

A  model  used  by  Parks  [Ref.  l:p.  4]  to  account  for  the 
rotating  atmosphere  separates  the  relative  velocity  of  the 
satellite  in  relation  to  atmosphere  into  two  parts.  The 
model  is  given  by: 
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Vrel  =  vt  /6   •  (11) 

Here  "v-^qi  is  the  satellite  to  atmosphere  relative  velocity, 
v^  is  the  satellite  tangential  velocity,  and  6  is  a  factor 
accounting  for  orbital  orientation  (inclination)  with 
respect  to  a  rotating  atmosphere  (see  Figure  11) . 

Equation  11  in  Parks  [Ref.  l:p.  4]  models  the 
inclination  factor  in  the  following  way: 

r'  2 

5   =   [1  -  (-^)  Aoj   cos  i]  (12) 

V     e 
P 

In  this  model,  6  is  the  rotation  factor,  rp  the  radius  at 
perigee,  Vp  the  satellite  velocity  at  perigee,  ooe  ^^^ 
earth's  angular  rotation  rate  (assumed  constant) ,  i  the 
orbital  inclination,  and  A  the  ratio  of  rates  of  rotation  of 
the  atmosphere  to  the  earth's  rotation. 

The  model  for  6  used  in  the  propellant  longevity  model 
neglects  the  effect  of  the  atmospheric  wind,  but  this  should 
introduce  only  a  small  error.  Recent  atmospheric  work, 
Jacchia  [Ref.  7],  suggests  a  significant  gravity  wave 
induced  motion  in  the  upper  atmosphere.  As  more  experimen- 
tal information  on  this  bulk  motion  becomes  available,  the 
incorporation  of  the  gravity  wave  wind  may  become 
increasingly  important. 

The  parameter  A  used  by  the  propellant  longevity  model 
in  the  satellite  velocity  multiplier  submodel  is  probably 
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not  constant.  Atmospheric  theory  suggests  it  may  be 
dependent  on  both  altitude  and  latitude.  The  range  of 
values  for  the  A  parameter  is  assumed  to  be  from  a  maximum 
of  about  1.0  (at  low  altitudes)  to  a  minimum  of  around  0.9 
(at  high  altitudes) . 

The  coding  of  the  theory  in  NSWC  TR  83-243  for  satellite 
fuel  mass  decrement  did  not  compute  6  internally.  It  was 
user  supplied  input  in  the  original  version.  As  part  of  the 
mass  decrement  model  validation  at  the  Naval  Postgraduate 
School,  code  internal  computation  of  5  was  added.  This 
allows  6  to  automatically  adjust  the  satellite-to-atmosphere 
relative  velocity  Vj-^i  for  the  results  predicted  by  the 
theory  as  perturbation-induced  changes  to  the  inclination 
occur.  For  a  long  mission  life  that  exhibits  significant 
inclination  drift,  the  internal  computation  of  should  add 
accuracy  to  the  mass  decrement  submodel.  In  place  of  the 
original  user  input  for  6  ,  the  A  is  now  specified  by  an 
input  file.   The  input  file  default  value  for  A  is  0.9. 

A  model  atmosphere  used  for  drag  computations  must  be 
able  to  correctly  predict  the  effects  on  density,  drag 
coefficient,  and  satellite-to-atmospheric  velocity.  There 
are  many  variations  in  the  atmosphere  that  affect  density. 
The  sum  of  all  these  effects  creates  a  density  bulge  in  the 
atmosphere,  the  location  and  shape  of  which  follow  maximum 
heating. 
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The  drag  coefficient  is  atmosphere  composition 
dependent.  Thus  the  model  atmosphere  must  be  able  to 
predict  the  dominant  species  in  the  satellite  environment. 

Satellite-to-atmosphere  relative  velocity  is  a  function 
of  satellite  tangential  speed,  atmospheric  rotation,  and 
orbital  inclination. 

The  Jacchia  J60  model  atmosphere  used  in  the  propellant 
longevity  model  only  includes  the  effects  of  diurnal  heating 
on  the  density  bulge.  There  are  no  computations  of 
atmospheric  composition.  The  newer  Jacchia  J77  model 
incorporates  all  the  known  density  variations  currently 
observed.  The  Jacchia  J77  atmosphere  also  predicts 
thermospheric  composition. 

The  background  concepts  for  the  propellant  longevity 
model  understanding  and  operation  have  been  presented.  The 
next  chapter  describes  the  model  and  its  computer  program. 
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III.   PROPELLANT  LONGEVITY  MODEL  DESCRIPTION 

This  chapter  describes  the  propellant  longevity  model. 
First,  a  brief  summary  of  the  model  theory  is  presented.  The 
second  part  of  this  chapter  describes  the  coding  of  the 
model  as  installed  at  the  Naval  Postgraduate  School. 
Finally,  the  third  section  gives  the  model  code  compiling 
and  handling  instructions  specific  to  the  system  at  the 
Naval  Postgraduate  School. 

A.   PROPELLANT  LONGEVITY  MODEL  THEORY 

A  detailed  development  of  the  propellant  longevity  model 
is  presented  by  Dr.  Parks  [Ref.  1].  The  propellant 
longevity  model  hinges  on  the  prediction  of  the  fuel  mass 
used  by  a  satellite  to  compensate  for  atmospheric  drag.  The 
mass  decrement  submodel  predicts  the  amount  of  fuel  consumed 
after  each  revolution  of  the  satellite. 

The  mass  decrement  submodel  for  a  low-earth-orbit 
satellite  in  an  oblate  diurnal  atmosphere  is  presented  in 
Reference  1  as  Equation  56.  The  diurnal  density  bulge  is 
the  only  density  and  compositional  atmospheric  variation 
considered  in  the  model. 

The  integrals  of  the  mass  decrement  model  (Equation  56) 
are  next  expressed  in  terms  of  Bessel  functions  of  the  first 
kind  and  order  n  as  defined  by: 
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27T 

1^(6  ae)   =   2^  /   cos(nE)  exp(3ae  oosE)  dE  (13) 


Here  I^  is  the  Bessel  function  of  the  first  kind  and  order 
n,  6  the  inverse  density  scale  height  (assumed  constant) ,  a 
the  semi-major  orbital  axis,  e  the  orbital  eccentricity,  n  = 
0,1,2,  ...,  6,  and  E  the  orbital  eccentric  anomaly.  All  of 
these  orbital  mechanics  terms  were  presented  in  Chapter 
II. A.   The  integral  is  to  be  evaluated  over  one  revolution. 

The  Bessel  function  expansion  of  the  mass  decrement 
model  is  given  in  equation  58  of  Reference  1.  This 
particular  expansion  is  central  to  the  model.  Provided 
accurate  values  for  other  multipliers  are  produced  by  the 
model  atmosphere,  (accurate  values  for  atmospheric  density, 
satellite  drag  coefficient,  and  satellite-to-atmosphere 
relative  velocity  are  essential)  greatly  increased  accura- 
cies should  result  in  predicting  the  fuel  consumed  to 
compensate  for  atmospheric  drag.  The  main  draw  back  to  the 
use  of  the  this  type  of  Bessel  function  is  the  known 
slowness  of  its  convergence. 

B.   COMPUTER  CODE  DESCRIPTION 

A  detailed  description  of  the  model  code  and  its  input 
elements  is  contained  in  the  code  and  included  in  Appendix 
A.   An  operational  description  is  presented  here. 

The  computer  code  calculates  the  propellant  mass 
remaining  together  with  the  cumulative  mission  life  after 


53 


each  revolution  of  the  satellite.   When  propellant  fuel  is 
exhausted,  the  computation  terminates. 

Computer  output  is  available  in  several  forms  and  is 
specified  by  the  user  through  the  use  of  some  control 
parameters.  The  user  can  select  one  of  two  modes;  standard 
mode  or  a  sensitivity  analysis  mode.  The  input  elements  used 
as  control  parameters  for  output  and  mode  control  are: 
maximum  number  of  revolutions  ITERl,  sensitivity  analysis 
granularity  PSIZE,  output  print  frequency  ITER2 ,  mode 
control  PSENS/  anomalistic  period  print  control  lET,  number 
of  scheduled  orbit  adjusts  NOA,  scheduled  orbit  adjust  array 
IRV(K)  ,  and  the  cunount  of  each  adjust  array  DA(K)  .  These 
terms  are  discussed  below. 

The  most  basic  form  of  output  available  is  a  tabulation 
of  propellant  remaining  and  accumulated  mission  life  for  20 
specified  satellite-orbital  initial  conditions.  The  20 
initial  conditions  are  specified  in  an  input  file. 
Accumulated  mission  life  is  expressed  in  modified  Julian 
days  (MJD) ^  and  revolution  number.  Frequency  of  output 
printing  is  selectable  from  a  maximum  of  once  per  revolution 
to  a  minimum  of  once  per  10000  revolutions  by  user 
specification  of  ITER2 . 

A  cap  on  runaway  computer  time  is  available  by  selection 
of  ITERl.   Satellites  with  small  cross-sectional  areas,  high 


^A  description  of  the  modified  Julian  day  time 
measurement  is  in  Taff  [Ref.  2:p.  103],  The  modified  Julian 
day  is  reckoned  from  0^  Nov.  17  1858. 
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altitudes,  and  large  amounts  of  maneuvering  fuel  can  have 
very  long  orbital  life-times.  By  defining  a  maximum  number 
of  revolutions  to  be  considered  in  the  calculations  through 
the  reasonable  selection  of  ITERl,  unnecessary  expenditure 
of  computer  time  can  be  avoided. 

As  part  of  the  model  verification  a  sensitivity  analysis 
option  was  added  to  the  program.  It  can  be  selected  by 
setting  PSENS  =  1.  The  standard  operating  mode  is  selected 
if  PSENS  =0.  In  the  sensitivity  analysis  mode,  the 
subroutine  SNALYS  creates  an  array  of  elements  selected  by 
the  user  to  be  linearly  scanned  from  a  specified  minimum  to 
maximum  with  a  given  granularity. 

Granularity  is  specified  by  the  selection  of  the  control 
parameter  PSIZE.  A  granularity  of  1/100  (code  specified  as 
PSIZE  =  100)  will  cause  100  executions  of  the  propellant 
longevity  code  to  complete  fuel  exhaustion.  This 
specification  can  lead  to  excessive  computer  execution  times 
for  high  altitude  satellites.  A  granularity  of  1/5  will 
cause  five  executions  of  the  code  to  fuel  exhaustion. 
However,  only  five  iterations  from  a  given  minimum  to  a 
given  maximum  may  be  too  coarse  to  reveal  the  true  trend  of 
the  satellite  mission  life. 

Each  time  the  code  is  sequentially  executed  in  the 
sensitivity  analysis  mode,  the  elements  selected  by  the  user 
for  analysis  are  incremented  linearly.  The  output  in  the 
sensitivity  analysis  mode  is  a  table  of  mission  lifetimes. 
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resulting  from  the  initial  values  of  input  elements  not 
selected  for  analysis,  and  the  current  value  of  the  elements 
selected  for  analysis. 

The  output  always  echoes  the  initial  input  element 
conditions  as  read  from  a  user  supplied  input  file.  Further 
choice  of  output  and  mode  control  parameters  allow  the  user 
to  tailor  the  output. 

Three  more  input  control  elements  allow  for  the 
specification  by  revolution  number  of  up  to  ten  adjusts  to 
the  semi-major  axis.  The  element  IRV(K)  is  an  array  whose 
elements  specify  the  orbital  revolution  number  for  axis 
adjust.  The  element  DA(K)  specifies  the  amount  of  change  in 
the  semi-major  axis  in  kilometers  for  each  scheduled 
adjustment.  Finally,  NOA  specifies  the  total  number  of 
orbital  adjusts  to  be  performed.  If  the  user  does  not 
desire  any  changes  to  the  semi-major  axis,  NOA  must  be  set 
to  0  and  no  values  for  IRV(K)  and  DA(K)  input  to  the  file 
stream. 

An  additional  twenty  input  elements  must  be  specified  by 
the  user  in  the  input  file.  Five  of  these  are  satellite 
physical  parameters,  five  are  atmosphere  specific 
parameters,  and  ten  are  the  Brouwer  mean  orbital  element  set 
of  the  NAVSPASUR  system. 

The  five  satellite  physical  parameters  are:  propulsion 
motor  specific  impulse  (SPI  in  seconds) ,  initial  mass  of  the 
satellite  including  fuel  (WT  in  kilograms) ,  initial  mass  of 
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the  maneuvering  fuel  (W  in  kilograms) ,  satellite  atmospheric 
drag  coefficient  (CD) ,  and  satellite  cross-sectional  area  (A 
in  square  kilometers) . 

Motor  specific  impulse  SPI  can  range  from  220  seconds 
for  monopropellants  to  3000  seconds  for  ion  drives.  The 
default  value  is  230. 

The  initial  satellite  mass  WT  depends  on  the  satellite 
mission  and  design.  A  typical  communications  satellite  has 
a  mass  of  about  3500  kilograms.  An  SDI  satellite  may  be 
considerably  more  massive.  The  default  value  for  the 
satellite  mass  is  9100  kilograms. 

Initial  maneuvering  fuel  mass  W  of  the  satellite  depends 
on  the  design  life  of  the  satellite  as  well  as  it's 
ballistic  coefficient.  Low  altitude  satellites  need  at 
least  one  percent  of  their  initial  mass  as  maneuvering  fuel 
if  mission  lives  of  years  are  to  be  achieved.  The  default 
value  for  initial  maneuvering  fuel  is  100  kilograms. 

The  drag  coefficient  CD  is  a  measure  of  the  slowing 
power  of  the  atmosphere  on  a  given  satellite.  In  the 
thermosphere  it  is  very  atmosphere  composition  dependent.  A 
value  of  2.2  has  been  settled  on  by  orbital  analysts  for 
altitudes  of  200  to  400  kilometers.  A  value  of  at  least  3.0 
must  be  assigned  at  altitudes  above  500  kilometers 
(corresponding  to  the  helium  atmosphere) .  The  default  value 
for  the  drag  coefficient  is  2.2. 
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Satellite  cross-sectional  area  is  again  very  design 
dependent.  The  echo  series  satellites  had  cross-sections 
measuring  in  hundred's  of  square  meters.  Some  SDI 
satellites  may  also  have  very  large  cross-sections.  The 
default  for  cross-sectional  area  is  50.0*10~^  square 
kilometers,  corresponding  to  50  square  meters. 

The  five  atmosphere  specific  parameters  are:  atmosphere- 
to-earth  relative  rotation  rate  D,  10.7  cm  solar  flux  F107 
described  in  Chapter  II,  ninety-day  average  solar  flux  FBAR, 
geomagnetic  activity  index  AKP,  and  the  time  of  vernal 
equinox  passage  TVE. 

The  atmosphere-to-earth  relative  rotation  rate  D  is  the 
factor  explained  at  the  end  of  Chapter  II.  The  default 
value  is  0.9. 

The  decimetric  solar  flux  and  the  average  decimetric 
solar  flux  F107  and  FBAR  are  again  explained  in  Chapter  II. 
Default  values  of  100  are  assigned  to  each. 

The  geomagnetic  index  AKP  is  discussed  in  Chapter  II. 
It  ranges  from  0  to  7  and  is  set  at  a  default  value  of  2.00. 

With  the  current  atmospheric  model  Jacchia  J60,  both 
FBAR  and  AKP  are  read  as  input  but  not  used  for  the  density 
computation.  As  seen  from  a  sensitivity  analysis,  these 
parameters  currently  have  no  effect  on  mission  life.  The 
Jacchia  J77  model  atmosphere  uses  these  now  stranded  input 
parameters.  The  future  incorporation  of  the  Jacchia  J77 
atmosphere  model  will  be  able  to  use  these  input  elements. 
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The  ten  NAVSPASUR  orbital  elements  are  the  final  input 
elements  in  the  standard  input  file.  They  are  all  explained 
in  the  orbital  mechanics  section,  and  in  the  code  included 
in  Appendix  A.  Briefly,  they  are: 

UJD:   time  of  epoch  in  modified  Julian  days 
ETU:   initial  mean  anomaly 
GO:   argument  of  perigee 

HO:   right  ascension  of  the  ascending  node 
ES:   eccentricity 
XI:   inclination 
AO:   semi-major  axis 
RND:   rate  of  change  of  mean  motion 
ESD:   eccentricity  decay  rate 
ADOT:   semi-major  axis  decay  rate. 
The  first  seven  of  these  are  orbital  elements  and  the  final 
three  are  the  results  of  perturbations.    The  NAVSPASUR 
orbital  element  set  is  the  result  of  orbital  measurements  of 
actual  satellites. 

The  final  block  of  input  in  the  input  file  is  the 
sensitivity  analysis  specification  array.  It  must  be 
included  only  when  the  sensitivity  analysis  mode  is 
specified.  The  array  consists  of  20  lines  used  to  make  up 
the  array  of  elements  to  be  scanned  by  the  program  when  the 
sensitivity  analysis  mode  is  selected. 

The  sensitivity  analysis  specification  array  is  divided 
into  four  columns.    The  first  column  contains  the  input 
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element  variable  name.  These  are  the  same  names  as 
specified  in  the  satellite  physical,  atmosphere  specific, 
and  orbital  elements  input  discussions. 

The  second  column  is  a  control  parameter.  When  set  to 
1.0,  the  corresponding  input  element  is  included  in  the 
sensitivity  analysis.  Columns  three  and  four  are  in  a 
Fortran  F16.9  format  and  they  specify  the  minimum  and 
maximum  range  of  the  element  to  be  scanned.  Defaults  for 
the  minimum  and  maximum  have  been  provided  in  specific  cases 
corresponding  to  logical  and  physical  ranges  of  the  input 
elements.  In  some  cases  zeros  appear  in  the  minimum  and 
maximum  columns.  These  zeros  correspond  to  input  elements 
not  physically  constrained.  An  example  of  an  input  element 
not  physically  constrained  is  the  time  of  epoch  (UJD) . 

If  a  1.0  is  entered  in  column  two  of  a  particular  row  of 
the  sensitivity  analysis  specification  array,  that  input 
element  is  included  in  the  analysis.  Any  number  of  the 
input  elements  may  be  selected  for  simultaneous  analysis. 
When  selecting  a  multivariable  sensitivity  analysis,  the 
user  must  understand  the  interactions  of  the  incrementing 
variables  and  ensure  that  a  solution  exists  in  the  mission 
life  range  being  considered. 

When  a  variable  is  included  in  the  sensitivity  analysis, 
the  incrementing  values  override  the  initial  conditions. 
All  input  elements  not  user  selected  for  analysis  (those 
with  zeros  in  the  second  column)  are  reset  to  their  initial 
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values  after  each  iteration  during  the  analysis.  Appendix  B 
contains  sample  input  element  files  and  their  corresponding 
outputs . 

Input  elements  and  output  forms  have  been  explained. 
The  program  flow  is  discussed  next.  The  16  subroutines  and 
eight  functions  comprising  the  propellant  longevity  model 
are  explained  in  Appendix  A  in  the  Fortran  code.  The 
routine  PLEP  initializes  the  system,  including  inputting 
program  control  parameters  and  initial  values  for  all 
elements.  If  the  sensitivity  analysis  mode  is  selected, 
subroutine  SNALYS  generates  the  array  to  be  scanned 
linearly,  complete  with  the  necessary  increments  for  each 
iteration.  Analysis  mode  causes  completion  of  all 
computations  to  propellant  exhaustion  PSIZE  times.  The 
routines  SNALYS  and  PLEP  are  completed  only  once  even  in 
analysis  mode. 

Program  control  is  transferred  from  the  initialization 
routine  PLEP  to  the  subroutine  KOZAI,  the  main  driver  and 
bookkeeper  of  the  system.  The  routine  KOZAI  first  causes 
the  NAVSPASUR  Brouwer  mean  element  set  to  be  c"onverted  into 
a  Keplerian  (Kozai-like)  mean  orbital  set,  by  calling  the 
subroutines  BRAUER,  PERIOD,  and  MEAN. 

Subroutine  MEAN  calls  functions  XSPA-XSPM  to  convert  the 
Brouwer  osculating  elements  generated  in  BRAUER  to  Keplerian 
orbital  elements.   The  Keplerian  elements  are  then  returned 


61 


to  the  subroutine  KOZAI.   In  standard  mode  the  mean  element 
set  is  given  as  output. 

All  of  the  above  computations  have  simply  conditioned 
the  input  elements  to  conform  to  the  requirements  of  the 
propellant  longevity  model.  The  routine  KOZAI  now  causes 
calculation  of  the  fuel  used  to  compensate  for  atmospheric 
drag  after  each  revolution.  The  computation  is  continued 
until  the  fuel  is  exhausted  or  the  maximum  number  of 
revolutions  (ITERl)  is  reached.  When  the  computation 
terminates,  program  control  is  transferred  back  to  PLEP. 

To  calculate  the  drag-compensating  fuel  required  for 
each  revolution,  KOZAI  calls  the  subroutines  PERIOD,  GEOP, 
THRST,  and  DRAG. 

The  subroutine  GEOP  calculates  the  change  in  orbital 
elements  due  to  the  aspherical  earth.  It  uses  the  J2 ,  J3 , 
and  J4  terms  in  a  Legendre  polynomial  expansion,  as 
explained  in  Chapter  II.  Aspherical  earth  calculations  are 
performed  once  each  revolution. 

The  subroutine  THRST  checks  each  revolution  for  a 
scheduled  orbital  adjust.  If  a  change  to  the  semi-major 
axis  is  scheduled,  it  calculates  the  change  to  the  axis,  the 
argument  of  perigee,  the  eccentricity,  and  the  mean  anomaly. 
It  also  computes  the  amount  of  maneuvering  fuel  required  to 
perform  the  adjustment.  All  changes  are  returned  to  the 
subroutine  KOZAI. 
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The  subroutine  DRAG  calculates  the  fuel  mass  needed  to 
overcome  the  atmospheric  drag.  This  calculation  is 
accomplished  by  calling  four  subroutines:  SOLOC,  SATLOC, 
SALT,  and  JAC60. 

The  routine  SOLOC  computes  the  right  ascension  and 
declination  of  the  sun  for  the  model  atmosphere  diurnal 
density  bulge.  The  routine  SATLOC  computes  the  right 
ascension,  declination  and  geomagnetic  latitude  of  the 
satellite.  This  is  accomplished  by  first  calling  the 
subroutine  POSVEL  to  calculate  satellite  position  and 
velocity.  The  routine  POSVEL  in  turn  calls  the  subroutine 
NWTRPH  to  solve  Kepler's  equation  (see  Chapter  II).  The 
routine  SALT  computes  the  satellite  altitude. 

Atmospheric  density  is  computed  by  the  subroutine  JAC60 
based  on  the  Jacchia  J60  model  atmosphere  discussed  in 
Chapter  II.  A  maximum,  minimum,  and  mean  density  that  the 
satellite  will  encounter  in  the  oblate  diurnal  atmosphere 
are  calculated. 

The  mass  decrement  equation  (Equation  58  of  Reference  1) 
is  contained  in  the  subroutine  DRAG.  The  function  MMBSIR  is 
used  to  perform  the  Bessel  function  expansion.  The  routine 
DRAG  returns  the  fuel  mass  decrement  due  to  atmospheric  drag 
computation  to  KOZAI .  The  routine  KOZAI  then  decrements  the 
remaining  fuel  by  this  amount  and  checks  for  fuel  exhaustion 
and  maximum  orbital  number.  If  fuel  is  less  then  zero,  or 
the  current  revolution  number  equals  ITERl,  calculations  are 
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terminated  and  control  is  returned  to  PLEP.    In  standard 
mode,  program  execution  now  halts. 

In  sensitivity  analysis  mode,  all  input  elements  not 
selected  for  analysis  are  returned  to  their  initial  values, 
elements  selected  for  analysis  are  linearly  incremented,  and 
KOZAI  is  again  called  by  PLEP.  The  subroutine  KOZAI  then  is 
executed,  as  previously  described,  with  the  new  values  for 
the  selected  elements.  Upon  fuel  exhaustion  or  reaching  the 
maximum  revolution  number,  KOZAI  terminates  calculation  and 
passes  program  control  back  to  PLEP.  This  process  continues 
until  the  elements  user  selected  for  analysis  have  been 
incremented  from  their  minimum  to  maximum.  When  the  maximum 
value  of  the  elements  user  selected  for  analysis  is  reached, 
program  execution  terminates  in  the  sensitivity  analysis 
mode . 

The  tabular  output  in  the  analysis  mode  is  printed  line- 
by-line  in  KOZAI  upon  reaching  fuel  exhaustion  or  maximum 
revolution  number.  The  result  is  a  table  of  data,  PSIZE 
long,  that  represents  the  linear  scan  of  the  specified 
elements.  The  first  line  is  the  mission  life  corresponding 
to  the  minimum  value  of  the  scanned  elements.  Each 
subsequent  line  represents  the  total  mission  life  resulting 
from  an  increment  of  the  elements  selected  for  analysis  and 
the  initial  conditions  of  the  elements  not  selected.  The 
final  line  is  the  mission  life  resulting  from  selected 
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elements  at  their  maximum  value  and  non-selected  elements  at 
their  initial  conditions. 

Mission  life  is  expressed  in  revolution  number  and 
modified  Julian  days.  The  last  column  of  the  tabular  output 
in  analysis  mode  represents  the  current  value  of  the 
incremented  input  element.  When  more  than  one  element  is 
analyzed,  the  last  column  is  the  current  fraction  of  scan,  a 
real  number  between  0  and  1. 

The  tabular  output  can  easily  be  converted  into  an  input 
file  for  a  graphics  routine  by  deleting  the  first  98  lines 
of  output. 

C.   MODEL  CODE  OPERATION 

This  section  describes  the  operation  of  the  computer 
code  at  the  Naval  Postgraduate  School  on  the  IBM  3  033 
computer  under  VMS  using  the  VS  FORTRAN  compiler. 

The  propellant  longevity  model  must  be  compiled  in  the 
VM  mode  using  the  automatic  precision  increase  option.  If 
the  propellant  longevity  code  is  in  filename  SDI  and 
filetype  FORTRAN  (the  filetype  must  be  FORTRAN) ,  the  correct 
form  for  the  compilation  command  is: 

FORTVS  SDI  (AUT0DBL(DBLPAD4) ) 
The  listing  generated  by  this  compilation  takes  up  about  2  0 
percent  of  the  user's  "A"  disk. 

Prior  to  program  execution,  the  input  and  output  files 
must  be  specified.  If  the  input  file  is  in  filename  PLEPOl 
and  filetype  DATA2,  then  the  correct  form  of  the  command  is: 
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FILEDEF  5  DISK  PLEPOl  DATA2  Al 
This  informs  the  computer  where  to  look  for  input  during 
execution  of  the  model  code.  If  hard  copy  of  the  output  is 
desired,  it  can  be  written  to  a  file  instead  of  the  screen 
by  using  the  FILEDEF  command.  If  the  file  name  is  SDI  and 
the  desired  type  is  DATA,  the  correct  command  is: 

FILEDEF  6  DISK  SDI  DATA  Al 
After  execution  of  the  propellant   longevity  code,   the 
results  will  be  in  SDI  DATA  on  the  user's  "A"  disk. 

To  execute  the  model  listing  generated  upon  compilation, 
the  command  RUN  SDI  is  typed  (this  assumes  that  SDI  is  the 
filename  of  the  propellant  longevity  program) . 

This  program  can  have  very  long  execution  times.  A 
general  rule  of  thumb  is  that  execution  time  takes  less  than 
PSIZE  *  ITERl/10,000  minutes.  Roughly,  the  computation 
takes  100  microseconds  for  each  revolution.  For  example,  if 
PSIZE  is  100  and  ITERl  is  10,000  (corresponding  to  a  five 
year  lifetime)  then  the  IBM  3033  will  take  200  minutes  to 
complete  the  calculations.  The  user  must  carefully  choose 
small  values  of  PSIZE  and  ITERl  on  the  initial  investigation 
of  elements. 

The  simplest  graphical  routine  to  construct  plots  of 
tabular  data  is  EASYPLOT.  The  EASYPLOT  routine  was  used  to 
generate  all  the  plots  shown  in  Chapter  IV. 

In  summary,  the  steps  of  program  operation  are: 

(1)   Compile  the  program  using  the  automatic 
precision  increase  option. 
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(2)  Select  the  mode  of  operation  and  input  element 
values  by  manipulation  of  the  input  file. 

(3)  Define  input  and  output  files. 

(4)  Execute  the  program  by  typing  the  file  name. 

The  next  chapter  presents  the  results  of  the  input 
element  sensitivity  analysis  perfoirmed  to  test  the  model 
reasonableness  during  verification. 
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IV.   MODEL  VERIFICATION 

This  chapter  presents  the  results  of  the  propellant 
longevity  model  verification  at  the  Naval  Postgraduate 
School.  Standard  model  verification  tests  the  model  against 
three  questions.  First,  does  the  model  address  the  problem? 
Second,  are  the  model  predictions  reasonable?  Third,  how  do 
the  model  predictions  compare  to  real-world  observations? 

This  verification  addresses  only  the  first  two  criteria. 
Observational  data  for  comparison  with  the  model  predictions 
is  not  available  in  an  unclassified  state.  The  quantitative 
analysis  of  model  predictions  is  thus  not  addressed  in  this 
verification.       ^ 

The  model  does  address  the  initial  problem.  It  predicts 
the  fuel  mass  decrement  of  a  low-earth-orbit  satellite  doing 
intrack  micro-thrusting  to  overcome  atmospheric  drag. 

The  reasonableness  test  is  the  subject  of  the  remainder 
of  this  chapter.  An  effective  way  to  test  the  model 
performance  is  to  examine  model  predictions  resulting  from 
the  full  "intended  use"  range  of  the  input  variables. 
Furthermore,  if  a  plot  is  made  of  predicted  results  as  a 
selected  input  element  is  incremented  through  its  range, 
prediction  trends  can  be  found  and  examined.  This  procedure 
provides  a  qualitative  analysis  of  model  sensitivity  to  the 
selected  input  element. 
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To  aid  in  these  tests  and  trend  studies,  a  sensitivity 
analysis  feature  was  incorporated  into  the  program  code  of 
the  model,  as  previously  described.  Graphical  results  of 
each  single  variable  sensitivity  analysis  are  presented  for 
the  input  elements. 

The  format  of  the  following  figures  are  all  related. 
They  represent  plots  of  satellite  mission  life  in  modified 
Julian  days  (here  after  referred  to  simply  as  "days")  versus 
a  specified  input  element.  All  figures  have  the  Y-axis  as 
mission  life  and  the  X-axis  as  the  variable  selected  for 
analysis.  In  each  case  the  range  of  values  for  the  input 
element  is  selected  to  remain  within  reasonable  values  of 
the  input  variable,  as  discussed  in  Chapter  III.  The 
granularity  in  most  cases  is  1/25  (PSIZE  =  25)  .  In  some 
cases  it  is  necessary  to  decrease  the  granularity  to  1/100 
(PSIZE  =  100)  to  see  the  actual  trends. 

The  default  values  discussed  in  Chapter  III  are  used  for 
all  input  elements  not  selected  for  analysis  (except  where 
noted) .  They  correspond  to  a  near-polar,  retrograde  orbit 
of  an  actual  satellite  supplied  as  sample  input  by  NSWC. 
The  model  predicted  propellant  longevity  using  all  default 
conditions  and  10  kilograms  of  fuel  is  13.74  days.  When  the 
default  values  and  100  kilograms  of  fuel  are  used,  the 
resulting  propellant  longevity  is  150.21  days.  An  *  symbol 
on  the  plotted  curves  indicates  the  default  value  of  the 
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analyzed  input  element.   The  default  value  names  are  again 
summarized  here.  The  units  can  be  found  in  Chapter  III. 
Satellite  physical  parameters: 

SPI  (motor  specific  impulse)   230 

WT  (satellite  initial  mass) 

W  (initial  fuel  mass) 

CD  (drag  coefficient) 

A  (cross-sectional  area) 
Atmospheric  specific  parameters: 

D  (atm.  rotation  factor)       1.0 

F107  (decimetric  solar  flux)   100 

FBAR  (average  F107) 

AKP  (geomagnetic  index) 

TVE  (time  of  vernal  equinox)   43222.7382 
NAVSPASUR  orbital  elements: 

UJD  (time  of  epoch) 

ETU  (mean  anomaly) 

HO  (R.A.  of  ascending  node) 

GO  (argument  of  perigee) 

ES  (eccentricity) 

XI  (inclination) 

RND  (time  rate  of  mean  motion)  .000783912 

ESD  (eccentricity  decay  rate)  -.239590000 

AO  (semi-major  axis)  1.06 

ADOT  (semi-major  axis  decay)   -701.8652 


100 
2.0 


44619.98716775 

150.0000 

95.0000 

200.0000 

.0100000 

95.0000 
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The  input  element  sensitivity  analyses  are  divided  into 
related  sets.  All  the  sets  are  discussed  in  the  same  basic 
format.  Each  analysis  is  single-variable  and  arranged  as 
follows: 

(1)  The   expected   trend   of   real-world   mission   life 
resulting  from  the  selected  variable 
incrementing  through  its  range. 

(2)  The  expected  mission  life  trend  predicted  by 
the  model  as  a  result  of  the  selected 
variable  incrementing  through  its  range. 

(3)  The  actual  mission  life  results  of  the 
selected  variable  incrementing  through  its 
range. 

(4)  A  discussion  of  the  selected  variable 
sensitivity  analysis. 

The  explanations  for  the  plotted  trends  are  based  on 
expected  behavior  in  agreement  with  current  theory.  An 
understanding  of  the  actual  mechanisms  producing  the 
observed  trends  will  be  improved  with  further  research. 

For  a  given  set  of  satellite-orbital  configurations,  the 
satellite  mission  life  will  be  shortest  for  satellites 
traveling  through  the  most  dense  atmosphere.  Thus  for  a 
given  input  variable  sensitivity  analysis,  the  range  of 
values  corresponding  to  satellite  travel  through  the  most 
dense  atmospheric  conditions  will  produce  the  shortest 
mission  life.  The  range  of  values  of  the  same  input 
variable  corresponding  to  satellite  travel  through  the  least 
dense  atmosphere  will  likewise  produce  the  longest  mission 
life.      The   orbital-perigee   atmospheric-density-bulge 
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encounter  geometry  is  thus  critical  in  mission  life 
analysis. 

The  initial  sensitivity  plot  of  an  input  variable  is 
accomplished  with  10  kilograms  of  fuel.  In  many  cases,  true 
geopotential-induced  orbital-drift  trends  are  not  seen  with 
small  amounts  of  drag-compensating  fuel.  To  verify  trends 
caused  by  geopotential-induced  orbital  drift,  the  initial 
compensating  fuel  is  increased  to  100  kilograms.  The 
differences  in  the  10  kilogram  and  100  kilogram  mission  life 
plot  trends  can  then  be  attributed  to  geopotential  effects 
and  the  movement  of  the  center  of  most  direct  solar  heating 
during  the  mission. 

The  133.14  day  nominal  mission  life  resulting  from  the 
default  values  of  the  input  elements  and  100  kilograms  of 
drag  compensating  fuel  is  almost  a  half-year.  During  this 
half-year,  the  earth  is  revolving  around  the  sun  causing  the 
noon-sun  atmospheric  position  to  increase  in  longitude  by 
about  one  degree  per  day.  This  greatly  complicates  trend 
analysis  when  using  the  100  kilogram  fuel  case  alone.  It  is 
thus  necessary  to  include  the  10  kilogram  fuel  case.  The 
13.74  nominal  mission  life  of  the  10  kilogram  fuel  case  is 
short  enough  so  that  it  shows  little  longitudinal  movement 
of  the  most-direct-sun  atmospheric  position. 

In  many  of  the  sensitivity  analysis  plots,  the  predicted 
mission  life  curve  has  a  "knee."  This  "knee"  corresponds 
to  the  range  of  the  selected  input  variable  that  produces 
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changing  behavior  in  the  mission  life  curve.  The  location 
of  these  ranges  has  significance  to  design  research. 

The  input  variables  are  divided  into  four  categories 
based  upon  their  effect  on  mission  life.  The  first  category 
contains  the  input  variables  that  affect  the  atmospheric- 
bulge  orbital-perigee  encounter  geometry.  The  second 
category  contains  the  input  variables  that  affect  the 
slowing  effectiveness  of  the  atmospheric  density  bulge  on 
the  particular  satellite  configuration.  The  third  category 
contains  the  input  variables  that  affect  the  depth  of 
penetration  of  the  satellite  into  the  atmosphere  at  perigee. 
Finally,  the  fourth  category  contains  the  input  variables 
that  have  no  affect  at  all  on  mission  life.  As  will  be 
seen,  predicted  mission  life  is  most  sensitive  to  the  third 
category  (factors  affecting  atmospheric  penetration  at 
perigee) . 

Factors  affecting  the  atmospheric-density-bulge  orbital- 
perigee  encounter  geometry  are:  the  time  of  vernal  equinox 
passage  TVE,  the  time  of  epoch  UJD,  the  argument  of  perigee 
GO,  the  right  ascension  of  the  ascending  node  HO,  and  the 
inclination  XI.  The  time  of  vernal  equinox  provides  a  real 
time  reference  for  earth  in  its  orbit  relative  to  the  sun. 
When  combined  with  the  time  of  epoch  (sometimes  called  the 
time  of  perigee  passage) ,  these  two  factors  determine  the 
position  of  the  noon  sun  in  right  ascension  and  declination 
relative  to  the  atmosphere.   The  default  values  used  in  the 


73 


model  verification  position  the  initial  noon  sun  at 
approximately  3  00  degrees  longitude  and  20  degrees  south 
latitude  relative  to  the  celestial  meridian  containing  the 
first  point  of  aries.  The  right  ascension  of  the  ascending 
node  and  the  argument  of  perigee  combine  to  position  the 
point  of  perigee  in  longitude  relative  to  the  first  point  of 
aries.  The  orbital  inclination  positions  the  point  of 
perigee  in  latitude  (or  declination) .  Actual  positioning  of 
the  satellite  in  the  above  orbit  is  accomplished  by- 
specifying  the  mean  anomaly.  All  these  orbital  mechanics 
terms  have  been  explained  in  Chapter  II. A. 

The  expected  real-world  behavior  of  satellite  mission 
life  as  a  result  of  varying  the  time  of  vernal  equinox  or 
time  of  epoch  is  discussed  next.  The  time  of  vernal  equinox 
TVE  serves  to  set  a  real-time  reference  for  positioning  of 
the  earth  in  its  orbit  relative  to  the  sun.  When  the  time 
of  epoch  UJD  is  specified,  the  position  of  the  point  of  the 
most  direct  atmospheric  heating  (noon  sun)  in  the  atmosphere 
is  set.  The  atmospheric  density-bulge  follows  the  location 
of  the  most  direct  solar  heating.  The  result  of  TVE  and  UJD 
varying  in  relation  to  each  other  is  a  seasonal  positioning 
of  the  earth  in  its  orbit  around  the  sun  at  the  time  of 
initial  perigee  passage.  Due  to  the  earth's  rotational-axis 
tilt,  the  point  of  most  direct  solar  heating  moves  in 
longitude  and  latitude  as  the  earth  moves  in  its  orbit 
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around  the  sun.   Figures  3  and  4  in  Chapter  II. A  aid  this 
visualization. 

The  earth  and  solar  cycles  effecting  the  atmospheric 
density-bulge  combine  to  define  its  size  and  shape.  These 
cycles  are  discussed  in  Chapter  II. B,  and  stated  again  here. 
They  are:  diurnal,  semi-annual,  seasonal-latitudinal,  solar 
rotational,  and  solar  activity  related.  The  Jacchia  J60 
model  atmosphere,  as  used  in  the  propellant  longevity  model, 
incorporates  only  the  diurnal  cycle. 

Once  the  position,  size,  and  shape  of  the  atmosphere  is 
set,  the  mission  life  of  a  drag  compensating  satellite 
orbiting  within  it  can  be  developed.  For  such  a  satellite 
in  a  specific  initial  orbit,  the  following  real-world  trend 
in  mission  life  (propellant  longevity)  should  be  observed. 
The  trend  should  be  cyclic,  resulting  from  the  superposition 
of  all  the  atmospheric  heating  cycles.  Mission  lives  will 
be  shortest  for  orbits  having  their  perigee  in  the  most 
dense  atmospheric  bulges.  The  longer  mission  lives  will  be 
observed  for  satellites  having  their  orbits  farthest  away 
from  the  density-bulge.  Thus  the  plots  of  satellite  mission 
life  will  show  peaks  and  valleys  consistent  with  the  above 
discussion. 

In  the  propellant  longevity  model,  the  diurnal  cycle  is 
the  only  atmospheric  heating  cycle  that  is  incorporated. 
Positioning  of  the  earth  in  its  orbit  is  accomplished  by  the 
combination  of  the  time  of  vernal  equinox  and  time  of  epoch 
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in  the  subroutine  SOLOC.  As  discussed  in  Chapter  III,  the 
subroutine  SOLOC  positions  the  sun  by  right  ascension  and 
declination  over  the  earth's  atmosphere.  The  routine  DRAG 
considers  the  diurnal  density  bulge  to  be  centered  under  the 
solar  position.  The  density  profile  of  the  atmospheric 
density-bulge  is  generated  in  the  routine  JAC60  based  on  the 
value  of  the  input  element  F107 . 

The  propellant  longevity  model  considers  the  atmosphere 
to  be  oblate,  diurnal,  and  revolving  around  the  sun  with  a 
23.5  degree  rotational-axis  tilt.  The  expected  model- 
predicted-mission-life  trend  as  is  cyclic  as  TVE  and  UJD  are 
varied.  Every  year  should  be  an  exact  replica  of  every 
other  year.  Two  peaks  and  two  valleys  in  the  yearly 
mission-life  plot  should  be  observed.  The  deepest  valley 
occurs  when  the  point  of  perigee  is  most  closely  centered  in 
the  atmospheric  density  bulge.  With  a  low  eccentricity  (e  = 
.01)  orbit  considered,  a  much  shallower  valley  will  occur 
about  six  months  later.  This  shallower  valley  corresponds 
to  the  point  of  apogee  being  most  closely  centered  in  the 
atmospheric  density-bulge. 

Two  peaks  in  plotted  mission  life  should  occur 
corresponding  to  the  orbital  pane  being  most  perpendicular 
to  the  earth-sun  line.  The  peaks  should  be  six  mgnths 
apart,  of  approximately  equal  magnitude,  and  interspaced 
between  the  mission-life  valleys.  As  the  eccentricity  of 
the  orbit  is  increased,  the  apogee-associated  valley  will 
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become  shallower  while  the  perigee-associated  valley  becomes 
deeper. 

Figures  12,  13,  14,  and  15  present  the  plots  of  the 
sensitivity  analysis  of  TVE.  In  Figure  12,  TVE  has  a  range 
of  1000  days,  from  38422  to  39422  days.  The  initial  fuel  is 
10  kg.  The  resultant  predicted  mission  life  ranges 
cyclically  between  11  and  18  days. 

Figure  13  has  an  incremented  range  for  TVE  of  38422  to 
38522,  or  100  days.  Initial  fuel  is  10  kg.  The  resulting 
predicted  mission  life  again  ranges  cyclically  between  11 
and  18  days. 

In  Figure  14  the  width  of  the  incremented  range  is 
raised  to  10,000  days  (about  27  years).  The  minimum  value 
of  TVE  is  38422  and  the  maximum  is  49422.  Again  10  kg  of 
fuel  is  used.  The  results  are  cyclic  between  11  and  18  days 
until  TVE  becomes  greater  than  time  of  epoch  UJD.  The 
results  at  the  high  end  of  TVE  in  this  plot  are  thus 
unreasonable.  The  time  of  vernal  equinox  passage  chosen 
must  occur  before  the  satellite  is  launched  requiring  TVE  to 
be  less  then  UJD. 

The  TVE  is  again  incremented  through  a  small  range  of  2  0 
days  in  Figure  15.  Initial  fuel  is  100  kg.  The  resulting 
range  of  mission  life  predicted  is  from  150  to  153  days. 
Note  the  absence  of  the  cyclic  trend". 

Figures  16  and  17  present  the  results  of  the  sensitivity 
analysis  plots  of  the  time  of  epoch  UJD  with  a  range  of  1000 
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Figure  12. 


Mission  Life  versus  Time  of  Vernal 
Equinox  (1000  Days) 


78 


oo. 


t^_ 


in 

s 

P3H 


o 

CO 
C/5 


o_ 


o_ 


TT  _ 


n 


M_ 


38420      38440  38460  38480  38500  38520 

TIME  OF  VERNAL  EQUINOX  (M.J.D.) 


Figure  13. 


Mission  Life  versus  Time  of 
Vernal  Equinox  (100  Days) 
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Figure  14 . 


Mission  Life  versus  Time  of 
Vernal  Equinox  (10,000  Days) 
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Figure  15. 


Mission  Life  versus  Time. of 
Vernal  Equinox  (20  Days) 
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Figure  16. 


Mission  Life  versus  Time  of  Epoch 
(1000  Days  and  10  kg  of  Fuel) 


82 


43223   43423      43623      43823      44023      44223 

TIME  OF  EPOCH  (MJD) 


Figure   17. 


Mission  Life  versus  Time  of  Epoch 
(100  Days  and  100  kg  of  Fuel) 
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days.  Note  the  similarity  to  the  TVE  plots,  as 
theoretically  expected.  Figure  16  shows  the  10  kg  fuel 
case,  and  Figure  17  the  100  kg  fuel  case. 

The  trends  observed  in  Figures  12  through  17  are  in 
agreement  with  expected  model  predictions.  The  plotted 
trends  may  have  significant  differences  from  real-world 
trends  due  to  the  limitations  of  the  Jacchia  J60  model 
atmosphere.  It  should  be  noted  that  the  default  values  of 
TVE  and  UJD  produce  a  solar  location  approximately  30 
degrees  longitude  past  the  winter  solstice  in  the  southern 
hemisphere  heading  toward  the  spring  (vernal)  equinox.  The 
approximate  latitude  of  the  default  solar  position  is  20 
degrees  south. 

The  right  ascension  of  the  ascending  node  HO  positions 
the  orbital  plane  with  reference  to  the  celestial  meridian 
containing  the  first  point  of  aries.  When  the  argument  of 
perigee  GO  and  the  inclination  XI  are  specified,  the  point 
of  perigee  is  positioned  in  longitude  and  latitude.  This 
geometry  is  shown  on  Figure  4  in  Chapter  II. A. 

The  real-world  trend  in  mission  life  of  a  drag 
compensating  satellite,  as  a  result  of  varying  HO,  GO,  and 
XI  is  cyclic.  The  cycle  repeats  every  3  60  degrees  for  HO 
and  GO,  and  every  180  degrees  for  XI.  The  observed  pattern 
of  peaks  and  valleys  will  be  highly  variable,  depending  on 
the  values  of  these  three  factors. 
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Model  predicted  mission  life  trends  due  to  varying  HO, 
GO,  and  XI  through  their  ranges  should  closely  parallel 
real-world  expectations.  This  close  parallel  between  model 
predictions  and  observation  assumes  the  atmospheric  density- 
bulge  has  been  properly  positioned  and  described. 

The  default  values  of  95  and  200  degrees  for  the 
inclination  and  the  argument  of  perigee  combine  to  produce 
a  point  of  perigee  that  is  very  nearly  20  degrees  below  the 
descending  node  in  the  southern  hemisphere  of  the  earth's 
atmosphere.  The  results  of  the  sensitivity  analysis  plots 
for  the  right  ascension  of  the  ascending  node  and  the 
argument  of  perigee  are  shown  in  Figures  18  through  21. 

The  right  ascension  of  the  ascending  node  HO  is 
incremented  from  0  to  3  60  degrees.  Figure  18  presents  the 
10  kg  fuel  case  and  Figure  19  the  100  kg  fuel  case. 
Predicted  mission  life  ranges  are  11  to  18  days  in  Figure 
18,  and  135  to  160  days  in  Figure  19. 

The  argument  of  perigee  GO  is  incremented  from  0  to  3  60 
degrees  with  a  10  kg  fuel  case  in  Figure  20,  and  a  100  kg 
fuel  case  in  Figure  21.  The  predicted  mission  life  are  13 
to  18  days,  and  148  to  160  days,  respectively. 

In  Figure  18  two  approximately  equal  mission  life  peaks 
and  valleys  are  observed.  The  valleys  correspond  to  the 
points  of  perigee  and  apogee  closely  approaching  the 
atmospheric  density  bulge.  The  peaks  occur  when  the  orbital 
plane  is  most  nearly  perpendicular  to  the  earth-sun  line 


85 


60  120  180  240  300 

RA  OF  ASCENDING  NODE  (DEC.) 


360 


Figure  18. 


Mission  Life  versus  Right  Ascension 
of  the  Ascending  Node  (10  kg  of  Fuel) 
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Figure  19. 


Mission  Life  versus  Right  Ascension 
of  the  Ascending  Node  (100  kg  of  Fuel) 
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Figure  20. 


Mission  Life  versus  Argument  of 
Perigee  (10  kg  of  Fuel) 
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Figure  21. 


Mission  Life  versus  Argument 
of  Perigee  (100  kg  of  Fuel) 
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which  contains  the  atmospheric  density  bulge  center.  The 
change  in  shape  of  the  10  and  100  kilogram  curves  (Figures 
18  and  19)  is  due  to  the  movement  of  the  atmospheric  density 
bulge  during  mission  life  and  to  geopotential-induced 
orbital  drift. 

With  default  values  of  95  degrees  for  both  the  right 
ascension  of  the  ascending  node  and  inclination,  the  effect 
of  varying  the  argument  of  perigee  from  0  to  3  60  degrees  is 
to  move  the  point  of  perigee  around  the  initially  fixed 
orbital  plane  (see  Figure  4  in  Chapter  II)  .  The  plot  of 
mission  life  corresponding  to  this  motion  should  by  cyclic 
showing  one  peak  and  one  valley.  The  valley  corresponds  to 
the  perigee  point  closest  to  the  atmospheric  density-bulge 
center.  The  peak  in  plotted  mission  life  occurs  when  the 
apogee  point  is  closest  to  the  atmospheric  density-bulge 
center. 

Figures  2  0  and  21  show  the  plotted  results  of  the 
sensitivity  analysis  of  the  argument  of  perigee  GO.  The 
expected  single  peak  and  single  valley  are  seen.  The 
difference  in  the  10  kilogram  case  in  Figure  20  and  the  100 
kilogram  case  in  Figure  21  is  due  to  the  seasonal  movement 
of  the  atmospheric  density-bulge  and  geopotential-induced 
orbital  drift.  The  magnitude  of  the  mission  life  variations 
is  less  than  that  produced  by  varying  the  right  ascension  of 
the  ascending  node.   This  effect  is  due  to  low  eccentricity 
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orbit  (e  =0.01)  and  the  position  of  the  orbital  plane  with 
reference  to  the  earth-sun  line. 

Factors  which  primarily  move  the  orbital  plane  in 
relation  to  the  atmospheric  density-bulge  center  will  show 
two  peaks  and  two  valleys  in  their  single-cycle  mission-life 
range.  Factors  which  primarily  move  the  point  of  perigee  in 
an  initially  specified  orbital  plane  will  show  one  peak  and 
one  valley  in  their  single-cycle  range.  These 
considerations  assume  that  only  low  eccentricity  orbits  are 
being  considered. 

Based  on  the  above  discussion,  Figures  18  through  21 
appear  to  agree  with  expected  real-world  and  model- 
predicted  mission  life  trends. 

The  analysis  of  predicted  mission  life  trends  due  to 
varying  the  inclination  XI  is  more  complex.  Several  factors 
must  be  considered  in  addition  to  the  orbital-perigee 
atmospheric  density-bulge  encounter  geometry.  First,  the 
diurnal  atmosphere  model  used  by  the  propellant  longevity 
model  is  also  oblate,  and  is  thus  more  dense  at  the  equator 
than  at  the  poles.  This  fact  means  that  satellites  in  near 
polar  orbits  travel  through  an  overall  less  dense  atmosphere 
than  those  in  a  near  equatorial  orbit.  Second,  the 
atmosphere  rotates  in  the  same  direction  as  the  earth. 
Hence  prograde  orbital  satellites  experience  less  satellite- 
to-atmosphere  .  relative  velocity  than  satellites  •  in 
retrograde  orbits.   If  the  satellite  altitude  is  low  enough. 
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a  near-equatorial  prograde  orbit  satellite  has  a  longer 
mission  life  than  the  same  satellite  in  a  near-equatorial 
retrograde  orbit.  These  additional  effects  complicate  the 
analysis  of  mission  life  resulting  from  varying  the  orbital 
inclination. 

The  default  value  position  of  the  center  of  the 
atmospheric  density-bulge  is  300  degrees  longitude  and  20 
degrees  south  latitude.  The  bulge  is  moving  at  a  rate  of 
about  one  degree  longitude  per  day  toward  the  vernal 
equinox.  The  point  of  perigee  is  default  positioned  at 
approximately  275  degrees  longitude  and  20  degrees  south 
latitude.  As  can  be  seen,  the  point  of  perigee  is  quite 
close  to  the  center  of  the  atmospheric  density-bulge.  The 
motion  of  the  point  of  perigee,  as  the  inclination  is 
increased  from  near  0  degrees  to  near  180  degrees,  describes 
a  semicircle  2  0  degrees  in  radius  around  the  descending 
node.  The  point  of  perigee  remains  quite  close  to  the 
initial  position  of  the  center  of  the  atmospheric  density- 
bulge. 

Figures  22,  23,  24,  and  25  are  the  predictions  of  the 
sensitivity  analysis  plots  of  orbital  inclination  XI.  The 
incremented  range  is  from  22.5  to  157.5  degrees.  Nearer 
equatorial  inclinations  are  prohibited  by  the  available 
arithmetic  precision  on  the  IBM  3033.  Figure  22  shows  the 
10  kg  fuel  case.  Figure  23  the  100  kg  fuel  case,  and  Figure 
24  the  100  kg  fuel  case  with  a  granularity  of  1/100  and  a 
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Figure  22. 


Mission  Life  versus  Orbital 
Inclination  (10  kg  of  Fuel) 
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Figure  23. 


Mission  Life  versus  Orbital 
Inclination  (100  kg  of  Fuel) 
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Figure  24. 


Mission  Life  versus  Orbital  Inclination 
(100  kg  of  Fuel  and  PSIZE  =  100) 
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Figure  25. 


Mission  Life  versus  Orbital  Inclination 
(100  kg  of  Fuel,  PSIZE  =  100,  and  D  =  0.9) 
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rotation  factor  of  0.9.  Predicted  mission  life  ranges  are 
7.8  to  13.7  days  in  Figure  22,  and  90  to  138  days  in  Figures 
23  and  25. 

Figure  25  shows  the  effect  of  the  using  100  kilograms  of 
compensating  fuel  and  a  rotation  factor  of  0.9.  This 
situation  is  compared  to  100  kilograms  of  compensating  fuel 
and  a  rotation  factor  of  1.0  in  Figure  23.  Only  very  slight 
differences  occur  between  Figures  23  and  25.  Mission  life 
is  a  little  longer  in  the  retrograde  orbits  using  a  rotation 
factor  of  0.9  and  a  little  shorter  in  the  prograde  orbits 
using  a  rotation  factor  of  0.9. 

When  the  satellite  orbit  is  lowered  to  an  altitude  at 
perigee  of  200  kilometers  (AO  =  1.03)  with  1000  kg  of  fuel, 
mission  life  is  reduced  in  the  retrograde-equatorial  orbit 
below  that  of  the  prograde-equatorial  orbit.  The  observed 
peaks  at  70  and  140  degrees  also  shift  leftward.  The 
observed  dimple  seen  in  Figure  22  at  the  critical 
inclination  of  63  degrees  may  be  geopotential-induced  or  an 
artifact  of  the  Kozai  approximations  as  discussed  by  Taff 
[Ref .  2:pp  332-342] . 

A  further  explanation  of  the  plot  of  mission  life  versus 
inclination  shown  in  Figure  22  is  presented  here.  A  more 
detailed  study  of  this  input  variable  is  .worthy  of  a 
complete  paper  itself. 

The  low  mission  life  found  at  low  values  of  inclination 
is  due  to  the  combination  of  the  point  of  perigee  position 
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and  oblate-atmosphere  effects.  The  low  inclination  orbits 
travel  through  the  "fat  part"  of  the  oblate  atmosphere  near 
the  equator.  As  the  inclination  increases,  the  point  of 
perigee  moves  away  from  the  center  of  the  atmospheric 
density-bulge  and  the  satellite's  orbit  moves  up  out  of  the 
thick  equatorial  atmosphere.  These  two  factors  combine  to 
produce  the  observed  rapid  rise  in  mission  life  seen  in 
Figure  22.  Beyond  near-polar  values  of  the  inclination,  the 
predicted  mission  life  again  falls.  The  mission  life  valley 
seen  at  120  degrees  is  probably  due  to  maximum  orbital  time 
near  the  atmospheric  density-bulge  center. 

Figure  23  shows  the  effect  of  starting  with  100 
kilograms  of  initial  drag-compensating  fuel.  Note  the  many 
peaks  and  valleys.  To  confirm  that  these  peaks  and  valleys 
represent  actual  predicted  mission  life  trends  and  not 
simply  data  scatter,  the  granularity  was  decreased  to  1/100. 
Figure  24  is  a  plot  of  these  results,  and  confirms  the 
trends  outlined  in  Figure  23.  The  many  peaks  and  valleys 
superimposed  over  the  general  trend  of  mission  life  seen  in 
Figure  22  are  probably  due  to  geopotential-induced  orbital 
drift. 

The  results  of  the  sensitivity  and  trend  analysis  of  the 
orbital  inclination  XI  have  been  discussed  with  reference  to 
in  Figures  22  through  25.  Future  investigation  of  the 
plotted  data  may  reveal  more  insight  toward  understanding 
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the  underlying  mechanisms.  The  sensitivity  analysis  plots 
show  no  trends  that  appear  unreasonable. 

The  sensitivity  analysis  for  the  input  variables 
effecting  the  orbital-perigee  atmospheric  density-bulge 
encounter  geometry  (category  one)  have  now  been  presented. 
The  next  set  of  analyses  are  those  of  the  second  category, 
which  contains  the  input  variables  affecting  the  slowing 
power  of  the  atmosphere.  These  variables  are:  decimetric 
solar  flux  F107,  drag  coefficient  CD,  cross-sectional  area 
A,  specific  impulse  SPI,  and  the  atmospheric  rotation  factor 
D.  The  initial  mass  of  drag-compensating  fuel  is  also 
included  in  this  category  because  it  improves  satellite 
mission-life  resistance  to  atmospheric  drag. 

The  real-world  expected  effect  on  mission  life  as  a 
result  of  varying  the  decimetric  solar  flux  F107  is  to  vary 
the  heating  effectiveness  of  the  sun.  As  the  decimetric 
solar  flux  increases,  the  density  of  the  atmosphere  will 
increase. 

The  Jacchia  J60  model  atmosphere  used  in  the  propellant 
longevity  model  subroutine  JAC60  uses  the  F107  flux  as  a 
multiplier  for  the  exponential  density  term.  The  expected 
model  predictions  of  mission  life  as  the  F107  flux  is  varied 
through  its  range  of  50  to  3  00  is  a  multiplicatively 
decreasing  plot. 

The  results  of  the  sensitivity  analysis  for  predicted 
mission  life  as  the  decimetric  solar  flux  is  varied  are 
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presented  in  Figures  2  6  and  27.  Figure  2  6  is  the  result  of 
10  kg  of  initial  maneuvering  fuel,  and  Figure  27  is  the 
result  with  100  kg  of  fuel.  The  predicted  mission  life  in 
Figure  26  ranges  from  24.99  to  3.96  days.  The  100  kg  fuel 
case  in  Figure  27  predicts  a  mission  life  ranging  from 
270.12  to  43.07  days.  These  results  are  in  agreement  with 
expectations . 

The  next  set  of  analyses  consider  effects  of  the  drag 
coefficient,  satellite  cross-sectional  area,  and  specific 
impulse.  All  of  these  factors  are  multipliers  acting  to 
change  the  resistance  of  the  satellite  to  atmospheric  drag. 
The  drag  coefficient  CD  is  a  proportionality  constant  in  the 
drag  force  model.  The  satellite  cross-sectional  area  A 
changes  the  size  of  the  satellite  that  the  atmosphere 
actually  acts  against.  As  satellite  cross-sectional  area 
increases  it  will  take  multiplicatively  more  energy  to 
compensate  for  the  orbital  energy  lost  by  the  satellite 
working  against  the  drag  force.  The  specific  impulse  SPI  is 
a  multiplier  that  expresses  the  effectiveness  of  the 
expelled  fuel  mass  in  overcoming  the  drag  force. 

The  rotation  factor  D  is  the  ratio  of  the  atmospheric 
rotation  rate  to  the  earth's  rotation  rate.  Its  effect  is 
to  change  the  satellite  to  atmospheric  relative  velocity  as 
the  inclination  changes.  A  logical  range  of  D  is  between 
the  values  0.9  and  1.0.  The  sensitivity  analysis  range  of 
0.1  to  10.0  was  chosen  to  show  model  sensitivity  for  a  wide 
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Mission  Life  versus  F10.7 
Flux  (10  kg  of  Fuel) 
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Figure  27. 


Mission  Life  versus  F10.7 
Flux  (100  kg  of  Fuel) 
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range  of  D.  As  presented  in  Chapter  II  Equation  12,  the  D 
factor  is  combined  with  the  inclination  and  the  ratio  of  the 
radius-at-perigee  to  the  velocity-at-perigee  producing  the 
6  factor.  This  factor  is  used  to  convert  the  satellite 
tangential  velocity  into  the  satellite-to-atmospheric 
relative  velocity. 

The  propellant  longevity  model  uses  A,  CD,  6  ,  and  SPI 
as  proportionality  constants,  multiplicatively  changing  the 
effectiveness  of  the  expelled  fuel  mass.  The  routine  PLEP 
combines  these  four  factors  into  a  single  constant  DEL  given 
by: 

DEL   =   A  *  CD  *  6   *  3.14159/SPI  (14) 

This  constant  DEL  is  passed  to  the  routine  KOZAI  which  in 
turn  passes  it  to  the  routine  DRAG.  In  the  routine  DRAG, 
DEL  is  used  as  the  final  multiplier  in  determining  required 
fuel  for  each  orbit. 

The  expected  model-predicted  mission-life  as  a  result  of 
changing  these  factors  is  a  multiplicative  increase  when  A, 
CD,  or  6  are  increased  linearly.  A  multiplicative  decrease 
should  occur  when  SPI  is  increased  linearly. 

Expected  results  of  varying  D  from  0.1  to  10.0  in  the 
near-polar  retrograde  orbit  are  shown  in  Figure  28.  Note 
that  a  very  slight  lowering  of  mission  life  is  seen  when  the 
rotation   factor   increases.   Figures  29   and  3  0   show  the 
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Figure  28. 


Mission  Life  versus  Rotation  Parameter 
(10  kg  of  Fuel  and  Default  Inclination) 
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Figure  29. 


Mission  Life  versus  Rotation  Parameter 
(10  kg  of  Fuel  and  Prograde  Orbit) 
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Figure  30. 


Mission  Life  versus  Rotation  Parameter 
(10  kg  of  Fuel  and  Retrograde  Orbit) 
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plotted  result  of  varying  D  through  the  same  range  of  0.1  to 
10.0.  However,  in  Figure  29  an  inclination  of  25  degrees  is 
used  instead  of  the  default  95  degrees,  and  in  Figure  3  0  an 
inclination  of  155  degrees  is  used.  In  Figure  29  it  is  seen 
that  the  mission  life  multiplicatively  increases  as  D  is 
linearly  increased.  Figure  30  shows  a  multiplicative 
decrease  in  mission  life  as  D  is  linearly  increased.  These 
trends  are  due  to  the  change  in  the  satellite-to-atmosphere 
relative  velocity.  The  increasing  mission  life  expectancy 
within  a  strongly  prograde  orbit  (Figure  29)  is  a  result  of 
the  decreasing  satellite-to-atmosphere  relative  velocity. 
The  decreasing  mission  life  expectancy  within  a  strongly 
retrograde  orbit  (Figure  30)  is  a  result  of  increasing 
satellite-to-atmosphere  relative  velocity.  The  plots  shown 
in  Figures  28,  29  and  30  are  in  agreement  with  expected 
results. 

The  results  of  the  sensitivity  analysis  plots  on  mission 
life  as  a  result  of  varying  CD,  A,  and  SPI  are  shown  in 
Figures  31  through  35. 

Figures  31  and  32  result  from  varying  the  satellite 
drag  coefficient  CD  from  1.0  to  3.5.  Figure  31  represents 
the  10  kg  fuel  case  and  Figure  32  the  100  kg  case.  The 
predicted  mission  life  range  in  Figure  31  is  31.18  to  8.63 
days,  and  in  Figure  32  from  3  08  to  99.58  days. 

The  sensitivity  analysis  mission-life  trend  plots  of 
satellite  cross-sectional  area  A  are  shown  in  Figures  3  3  and 
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Mission  Life  versus  Drag 
Coefficient  (10  kg  of  Fuel) 
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Figure  32. 


Mission  Life  versus  Drag 
Coefficient  (100  kg  of  Fuel) 
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Mission  Life  versus  Cross-Sectional 
Area  (10  kg  of  Fuel) 
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Figure  34. 


Mission  Life  versus  Cross-Sectional 
Area  (100  kg  of  Fuel) 
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Figure  35. 


Mission  Life  versus  Specific 
Impulse  (10  kg  of  Fuel) 
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34.  Figure  33  is  the  10  kg  fuel  case  and  Figure  34  the  100 
kg  case.  In  both  cases,  A  was  varied  from  0.000010  to 
0.0001000  square  kilometers.  The  ranges  mapped  by  these 
plots  are  80  to  0.64  days  and  736.29  to  108.7  days, 
respectively. 

Figure  35  is  the  analysis  of  motor  specific  impulse  SPI 
varying  from  200  to  3000  with  10  kg  of  initial  fuel.  The 
apparent  linear  shape  of  the  SPI  plot  is  due  to  its  inverse 
multiplicative  action  and  the  range  of  the  constant  DEL  (14) 
over  which  it  acts.  As  SPI  increases,  the  drag  compensating 
fuel  is  more  effective.  Linearly  increasing  SPI  has  the 
same  effect  as  starting  at  the  default  value  *  on  the  cross- 
sectional  plots  and  moving  along  the  plot  to  the  left.  Note 
that  the  shape  of  cross-sectional  mission  life  plot  in  this 
range  is  closer  to  the  expected  mirror  image  of  the  SPI 
plot. 

As  noted  earlier,  mission  life  plots  showing  "knees" 
have  significance  to  research  and  design  studies.  The  knees 
in  the  cross-sectional  area  mission  life  plots  are  very 
pronounced.  Less  distinctive  knees  are  seen  in  the  mission- 
life  plots  of  the  decimetric  solar  flux  and  the  drag 
coefficient.  The  default  value  of  the  cross-sectional  area 
is  right  in  the  knee  range  of  the  predicted  mission-life 
plots,  as  seen  in  Figures  34  and  35. 

The  final  input  variable  in  the  second  category  is  the 
initial  propellant  mass  W.   The  expected  mission  life  of  a 
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drag  compensating  satellite  as  fuel  is  increased  linearly 
should  be  a  linear  increase  in  the  predicted  mission  life. 
A  given  amount  of  fuel  produces  a  certain  mission  life. 
Twice  that  amount  of  fuel  in  the  same  general  configuration 
should  produce  about  twice  the  mission  life.  The  expected 
shape  of  the  sensitivity  analysis  plot  of  predicted  mission 
life  for  linearly  increasing  initial  propellant  mass  is  a 
approximately  straight  line.  Figure  36  is  the  result  of 
varying  W  on  mission  life.  The  plot  is  approximates  a 
straight  line  as  expected.  The  ripples  are  probably  due  to 
geopotential-induced  orbital-drift  or  to  movement  of  the 
atmospheric  density-bulge  during  mission  life. 

The  results  of  the  sensitivity  and  trend  analysis 
mission-life  plots  for  CD,  A,  D,  SPI,  and  W  agree  with 
expectations.  All  input  variables  in  the  second  category 
have  been  examined.  The  next  input  variables  to  be  examined 
are  in  the  third  category. 

Consider  now  the  input  variables  that  affecting  the 
depth  of  penetration  of  the  satellite  into  the  atmosphere  at 
perigee.  The  variables  that  affect  the  satellite  altitude 
at  perigee  are  the  semi-major  axis  and  the  orbital 
eccentricity.  These  variables  show  the  most  pronounced 
effects  on  mission  life.  The  real-world  effect  on  mission 
life  when  linearly  changing  the  altitude  at  perigee  is  an 
exponential  increase  in  mission  life.   The  effect  of  varying 
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Mission  Life  versus 
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the  eccentricity  and  the  semi-major  axis  on  the  altitude  at 
perigee  is  summarized  by 

rp  =  a  *  (1  -  e)  (15) 

Here,  rp  is  the  radius  at  perigee,  a  the  semi-major  axis, 
and  e  the  eccentricity.  A  linear  increase  in  the  semi- 
major  axis  results  in  a  modified  exponential  increase  in 
mission  life.  A  linear  increase  of  the  eccentricity  results 
in  a  modified  exponential  decrease  in  mission  life. 

The  propellant  longevity  model  computes  the  satellite 
altitude  at  perigee  in  the  subroutine  SALT.  The  subroutine 
DRAG  then  passes  the  altitude  at  perigee  to  the  routine 
JAC60  where  it  is  used  to  develop  an  exponent  for  an 
empirical  model  of  atmospheric  density. 

Mission  life  results  of  varying  the  orbital  eccentricity 
ES  from  0.01  to  0.1  are  presented  in  Figures  37  and  38. 
Figure  37  is  the  10  kg  fuel  case,  and  Figure  38  the  100  kg 
fuel  case.  The  predicted  mission  life  ranges  from  14  to  0 
days  in  Figure  37,  and  150  to  0  days  in  Figure  38. 

An  analysis  of  the  semi-major  axis  AO  is  presented  in 
Figure  39  with  a  range  of  1.025  to  1.175  earth  radii.  In 
this  case,  1  kg  of  fuel  is  allocated  initially.  The 
resulting  mission  life  is  from  0  to  384  days. 

The  plots  of  mission  life  resulting  from  varying  the 
semi-major  axis  AO  and  the  orbital   eccentricity  ES,   show 
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major  knees  in  their  curves.  The  *  symbol  indicates  the 
default  value  of  the  varied  element.  Note  that  the  amount 
of  fuel  used  in  the  analysis  of  the  semi-major  axis  is  1  kg. 
The  use  of  the  usual  10  kg  produces  a  plot  that  exceeding 
the  reasonable  range  of  mission  life.  The  results  shown  in 
Figures  37  through  39  are  in  agreement  with  expectations. 

The  input  variables  in  the  fourth  category  (no  effect  on 
mission  life)  are  the  mean  anomaly  ETU,  initial  satellite 
weight  WT,  average  decimetric  solar  flux  FBAR,  and  the 
geomagnetic  index  AKP.  FBAR  and  AKP  have  no  effect  on 
mission  life  due  to  the  limitations  of  the  Jacchia  J60 
atmosphere  model.  In  the  real  world  these  two  factors  may 
combine  with  the  decimetric  solar  flux  F107  to  produce 
changes  in  the  atmospheric  density  of  three  orders  of 
magnitude. 

The  model  theory  used  by  Dr.  Parks  [Ref.  1]  to  predict 
the  propellant  longevity  of  a  satellite  doing  intrack  micro- 
thrusting  to  overcome  atmospheric  drag  is  independent  of  the 
satellite  mass.  The  fuel  required  to  overcome  drag-induced 
orbital  energy  loss  is  based  strictly  on  the  satellite's 
shape  and  size. 

Figure  4  0  presents  the  sensitivity  analysis  plot  of 
predicted  mission  life  as  a  result  of  varying  the  satellite 
initial  mass  WT.  The  plot  is  a  straight  horizontal  line,  as 
expected. 
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The  final  analysis  is  concerned  with  the  mean  anomaly. 
In  the  near-circular  default  orbit,  the  mean  anomaly  and  the 
eccentric  anomaly  are  nearly  the  same.  The  Keplerian  mean 
anomaly  functions  to  position  the  satellite  in  its 
designated  orbit.  For  a  mission  life  lasting  several  days, 
it  should  make  no  difference  where  the  satellite  is 
initially  positioned  in  its  orbit.  Rather,  it's  the  orbit 
itself  that  makes  the  difference. 

Analysis  plots  of  the  m.ean  anomaly  ETU  are  presented  in 
Figures  41  and  42.  The  incremented  range  is  from  0  to  179 
degrees.  Figure  41  is  the  10  kg  fuel  case  and  Figure  42  is 
the  100  kg  fuel  case.  Predicted  results  range  from  13.74  to 
13.87  days  in  Figure  41,  and  from  149  to  154  days  in  Figure 
42.  ^ 

The  results  in  Figure  41  are  as  expected,  an  approximate 
horizontal  line.  In  Figure  42,  a  valley  is  seen  in  the  40 
to  120  degree  range.  The  valley  is  only  four  days  deep,  but 
is  difficult  to  explain.  Further  research  may  reveal  its 
underlying  mechanism. 

The  results  of  the  trend  and  sensitivity  analysis  have 
been  presented.  Mission  life  has  been  found  to  be  most 
sensitive  to  those  factors  that  affecting  the  orbital  radius 
at  perigee.  The  main  unexpected  results  are  found  in  the 
analysis  of  the  inclination  and  the  mean  anomaly.  The 
pattern  of  many  peaks  and  valleys  found  in  the  100  kilogram 
fuel  case  (Figures  23,  24  and  25)  warrants  further  research. 
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The  four  day  dimple  seen  in  the  mean  anomaly  plot  (Figure 
41)  also  warrants  further  study. 
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V.   SUMMARY 

This  final  chapter  presents  the  summary  and  conclusions 
of  the  thesis.  Recommendations  for  model  use  are  also 
given. 

The  low-earth-orbit  satellite  propellant  longevity  model 
has  been  tested  in  a  sensitivity  analysis.  Based  on  the 
results  of  the  sensitivity  analysis  the  model  reasonableness 
is  confirmed.  The  first  chapter  presented  a  model  overview, 
model  motivation,  and  the  process  employed  to  verify  the 
model.  Chapter  I  also  discusses  some  changes  made  to  the 
model  computer  program  for  the  IBM  3033  at  the  Naval 
Postgraduate  School.  Chapter  II  presents  the  orbital 
mechanics  and  atmospheric  concepts  as  background  for  model 
understanding  and  it's  operation.  In  Chapter  III  the  model 
and  its  computer  program  were  discussed.  Operating 
instructions  specific  for  the  computer  program  use  at  the 
Naval  Postgraduate  School  are  also  presented  there.  The 
sensitivity  analysis  results  are  presented  in  Chapter  IV. 

The  propellant  longevity  model  predicts  how  long  the 
fuel  of  a  low-earth-orbital  satellite  thrusting  to  overcome 
atmospheric  drag  will  last  before  starting  uncontrolled 
decay  into  the  atmosphere.  Model  verification  reasonable- 
ness tests  reveal  no  inconsistencies  with  the  established 
theory.     The   results   of  the  sensitivity  analysis   are 
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therefore  in  agreement  with  the  expectations.  The 
irregularities  in  the  model  predictions  (occurring  with 
inclinations  of  less  than  22.5  degrees)  are  due  to  computer 
round  off  error,  not  the  model  itself. 

The  propellant  longevity  model  currently  uses  the 
Jacchia  J60  model  atmosphere.  This  atmosphere  model  may 
significantly  reduce  the  real-world  accuracy  of  predicted 
mission  lifetimes. 

The  short  term  confidence  factor  in  predicted  mission 
life  might  be  increased  by  the  employment  of  the  Jacchia  J77 
model  atmosphere.  One  drawback  to  the  use  of  the  Jacchia 
J77  model  atmosphere  is  that  it  may  double  the  computer  run 
time. 

One  solution  to  the  trade-off  of  increased  accuracy  of 
the  Jacchia  J77  model  versus  the  shorter  computation  time  of 
the  Jacchia  J60  model,  would  be  to  allow  the  user  to  select 
an  atmosphere  option. 

Transportation  of  the  computer  program  of  the  model  to 
an  IBM  AT  could  provide  more  flexibility  in  preliminary 
design  analysis. 

This  model  is  a  tool  for  the  low-earth-orbit  satellite 
system  planner  in  designing  and  managing  satellite  systems 
that  have  precise  orbital  element  control  requirements. 
These   systems   are   those   with   SDI   and   reconnaissance 
missions. 
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This  model,  with  the  addition  of  the  Jacchia  J77  model 
atmosphere  and  a  subroutine  to  calculate  fuel  required  to 
compensate  for  aspherical  geopotential  forces,  could  be  used 
to  optimize  SDI  defensive  satellite  constellation  design. 

Existing  real-world  data  to  compare  with  the  model 
predictions  is  limited  to  a  very  classified  system.  The 
final  step  in  verifying  the  propellant  longevity  model  would 
be  to  make  this  comparison  test. 
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APPENDIX  A 
COMPUTER  PROGRAM  LISTING 


XX  XX 

XX  PROPELLENT  LONGEVITY                       XX 

XX  FOR                                 XX 

XX  LOW-ALTITUDE  EARTH  ORBIT  SATELLITES               XX 

XX  (A.D.  PARKS)                           xx 

X5<  XX 

X 
X 
X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


X 

X        THIS  PROGRAM  IS  THE  CODING  OF  THE  THEORY  CONTAINED  IN  NSHC 
X    TR  83-243.  IT  NAS  USED  TO  GENERATE  THE  NUMERICAL  EXAMPLES  ON 
X    PAGES  13  TO  17.    INITIALLY  CODED  FOR  THE  CDC  6700  IN  FORTRAN 
X    IV  BY  A.D.  PARKS,  THE  CODE  HAS  BEEN  TRANSPORTED  TO  THE  IBM  3033 
X    AND  TRANSLATED  TO  FORTRAN  77  BY  LT.  CD.  NOBLE.   ADDITIONAL  DOC- 
X    UMENTATION  WAS  ADDED  DURING  THE  TRANSPORTATION  AND  TRANSLATION. 
X 

X       THE  REFERENCE  MATERIAL  CITED  IN  THE  CODE  IS  MOSTLY  FORM  THE 
X    NAVAL  SURFACE  WEAPONS  CENTER  (NSWC)  TECHNICAL  REPORTS  (TR) . 
X    NSWC  TR  »  ARE  CITED  IN  THE  CODE  AS  COMMENTS.  OTHER  REFERENCES 
X    ARE  SITED  BY  NUMBER.  THEY  ARE: 
X 

REFERENCE  1       BROUWER,  D.,  "SOLUTION  OF  THE  PROBLEM  OF 

ARTIFICIAL  SATELLITE  THEORY  WITHOUT  DRAG'S 
THE  ASTRONOMICAL  JOURNAL,  VOL  64,  NO  1274, 
PP.  378-397,  1959. 


REFERENCE  2       AGRAWAL,  B.N.,  "DESIGN  OF  GEOSYNCHRONOUS 
SPACECRAFT",  PRENTICE-HALL,  1986. 

REFERENCE  3       JACCHIA,  L.G.,  "THERMOSPHERIC  TEMPERATURE, 
DENSITY  AND  COMPOSITION:  NEW  MODELS", 
SPECIAL  REPORT  375,  SMITHSONIAN  ASTRO- 
PHYSICAL  OBSERVATORY,  1977. 

REFERENCE  4       JACCHIA,  L.G.,  "EMPIRICAL  MODELS  OF  THE 
THERMOSPHERE  AND  REQUIREMENTS  FOR  IMP- 
ROVEMENTS", AFGL-TR-81-0195,  ADA101946, 
1981. 

REFERENCE  5       JACCHIA,  L.G.,  "ANALYSIS  OF  DATA  FOR  THE 
DEVELOPMENT  OF  DENSITY  AND  COMPOSITION 
MODELS  OF  THE  UPPER  ATMOSPHERE",  AFGL-TR- 
81-0230,  AD106420,  1981. 

REFERENCE  6       ORR,  L.H.,  "ORBITAL  LIFETIME  PROGRAM", 
NASA  TM  87587,  19S5. 

REFERENCE  7       TAFF,  L.G.,  "CELESTIAL  MECHANICS:  A  COMP- 
UTATIONAL GUIDE  FOR  THE  PRACTIOHER", 
JOHN  WILEY  S  SONS.  INC.,  1985. 

REFERENCE  8       FLEAGLE,  R.G.,  "AN  INTRODUCTION  TO  ATMOS- 
PHERIC PHYSICS",  ACADEMIC  PRESS,  1980. 
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MODULE 


PLEP 


SNALYS 


KOZAI 


BRAUER 


POSVEL 


NWTRPH 


FUNCTIONS 
XSPA-XSPR 


MEAN 


GEOP 


FUNCTION 


SYSTEM  INITIALIZATION  AND  CONTROL 
CALLED  BY:  CALLS: 

(SUBROUTINES) 
KXXX  SNALYS 

KOZAI 


TIMES  CALLED 
1 
1/ITERATION 


SENSITIVITY  MODE  SCANNED  ELEMENT  GENERATOR. 
CALLED  BY:  CALLS: 

(SUBROUTINES)   TIMES  CALLED 
PLEP  xxxx        X3(KX 

SYSTEM  WORKHORSE,  DRIVER  AND  RECORD  KEEPER. 
NSWC  TR  83-135. 
CALLED  BY: 

PLEP 


CALLS: 
(SUBROUTINES) 

BRAUER 

PERIOD 

MEAN 

THRST 

GEOP 

DRAG 


TIMES  CALLED 
1 
1 
1 

1/(SAT.REV, 
1/(SAT.REV, 


1/(SAT.REV.) 


CONVERTS  NAVSPASUR  BROUWER  MEAN  ORBITAL  ELEMENT 
SET  INTO  A  BROUWER  OSCULATING  ORBITAL  SET. 
A  CODING  OF  THE  BROUWER-LYDDANE  THEORY. 
NSWC  TR  83-107,  83-155,  BROUWER  (REF.  1). 


CALLED  BY. 
KOZAI 


CALLS: 
(SUBROUTINES) 


TIMES  CALLED 


SOLOC 


COMPUTES  SATELLITES  POSITION  AND  VELOCITY  IN 
INERTIAL  SPACE  FROM  THE  OSCULATING  ORBITAL 
ELEMENTS.  NSWC  TR  82-387  PP.  24-25. 
CALLED  BY:  CALLS: 

(SUBROUTINES)   TIMES  CALLED 
SATLOC  NWTRPH       1 

A  NEWTON-RAPHSON  SOUUTION  TO  KEPLERS  EQUATION. 
NSWC  TR  82  387  P  23. 
CALLED  BY:  CALLS: 

(SUBROUTINES)   TIMES  CALLED 
POSVEL  xxxx         xxxx 

COMPUTES  SHORT  PERIODIC  EFFECTS  ON  ORBITAL 
ELEMENTS.  CALLED  BY  SUBROUTINE  MEAN.  NSWC  81-456 
PAGES  9-11  AND  NSWC  TR  83  243  REF  3. 
CALLED  BY:  CALLS: 

(SUBROUTINES)   TIMES  CALLED 
MEAN  xxxx         xxxx 

USES  THE  WALTER  METHOD  TO  CONVERT  OSCULATING 
ORBITAL  ELEMENTS  TO  KOZAI  LIKE  MEAN  ELEMENTS 
USING  FIRST  ORDER  PERIODIC  VARIATIONS. 
SIMILAR  TO  THE  METHOD  IN  MSWC  TR  83  135  PP  16-18 
CALLED  BY:  CALLS: 

(FUNCTIONS)     TIMES  CALLED 
KOZAI  XSPA-XSPR     1 

COMPUTES  AVERAGED  SECULAR  CHANGES  IN  KOZAI  MEAN 
ELEMENTS  INDUCED  BY  GRAVITATIONAL  FIELD  THROUGH 
FOURTH  ZONAL  HARMONIC  (J2-J4  TERMS).  NSWC  TR  81- 
456.  EQS.(2.5). 
CALLED  BY:  CALLS: 

(SUBROUTINES)   TIMES  CALLED 
KOZAI  XKXX         X5(x» 

COMPUTES  THE  RIGHT  ASCENSION  AND  THE  DECLINATION 
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PERIOD 


DRAG 


JAC60 


THRST 


SALT 


SATLOC 


FUNCTION 
(MMBSIR) 


OF  THE  SUN  FOR  THE  MODEL  ATMOSPHERE 

456  PP.  20-21. 

CALLED  3Y:  CALLS: 

(SUBROUTINES) 
DRAG  XX3<?« 

COMPUTES  THE  ANOMALISTIC  PERIOD. 
CALLED  BY:  CALLS: 

(SUBROUTINES) 
KOZAI  xxxx 


NSWC  TR  81- 


TIMES  CALLED 


TIMES  CALLED 


COMPUTES  DRAG  EFFECTS  AND  FUEL  MASS  DECREMENT 
REQUIRED  TO  MAINTAIN  THE  SEMIMAJOR  AXIS  VIA 
IN  TRACK  MICROTHRUSTING.  NSWC  TR  82-2«+3  EQ . 
58.  THIS  ONE  EQUATION  IS  THE  HEART  OF  THE  MODEL 
CALLED  BY:  CALLS: 

(SUBROUTINES) 
KOZAI  SOLOC 

SATLOC 
SALT 
JAC60 
(FUNCTIONS) 

MMBSIR        1 
XSPM  1 


TIMES  CALLED 
1 
3 
1 
3 


COMPUTES  ATMOSPHERIC  DENSITY  BASED  ON  THE  JACCHIA 
(60)  MODEL  ATMOSPHERE.  CALLED  THREE  TIMES  TO 
COMPUTE  DIFFERENT  VALUES  FOR  THE  DENSITY  FOR  EACH 
REVOLUTION:  MAX,  MIN  AND  OP. 
CALLED  BY:  CALLS: 

(SUBROUTINES)   TIMES  CALLED 
DRAG  xxxK  xxxx 

CALCULATES  CHANGES  TO  KOZAI  ORBITAL  ELEMENTS  IF 
AN  ORBITAL  ADJUST  IS  SCHEDULED.  CALCULATES  THE 
MASS  OF  FUEL  USED  TO  PERFORM  THE  ADJUST. 
NSHC  TR  83-31  EQ .  16,  AND  NSWC  TR  81-456  EQ .  2-26 
CALLED  BY:  CALLS: 

(SUBROUTINES) 
KOZAI  xxxx 


TIMES  CALLED 


COMPUTES  SATELLITE  ALTITUDE.  NSWC  TR  81-456  EQ 

4.34. 

CALLED  BY:  CALLS: 

(SUBROUTINES) 
DRAG  xxxx 


TIMES  CALLED 


COMPUTES  SATELLITE  LOCATION.  RIGHT  ASCENSION, 
DECLINATION  AND  ITS  GEOMAGNETIC  LATITUDE. 
CALLED  BY:  CALLS: 

(SUBROUTINES) 
DRAG  POSVEL 


TIMES  CALLED 
1 


COMPUTES  BESSEL  FUNCTIONS.  AN  IM3L  ROUTINE 
RESIDENT  ON  THE  IBM  3800.  THIS  WILL  VARY  WITH 
LOCAL  INSTALLATIONS. 
CALLED  BY:  CALLS: 

(SUBROUTINES)   TIMES  CALLED 
DRAG  XX5(X  xxxx 


XXXXXXXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXKXXXXXXXXXXXXXXXXXXXXXXXKXXXXXXXX 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

XX  XX 

XX       PROGRAM  PLEP         XX 
XX  XX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
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PROGRAM  PLEP 
X 

X  THIS  MODULE  IS  THE  INPUT,  CONTROL  AND  INITIALIZATION  SECTION  FOR  THEx 

X  REST  OF  THE  SYSTEM.  IT  ACCEPTS  THE  'FIVE  CARD'  FORMAT  OF  MAVSPASUR   X 

K  AS  THE  BROUWER  MEAN  ORIBITAL  ELEMENT  SET.   THE  SPECIFIC  SATELLITE    x 

n  PARAMETERS  ARE  ENTERED.  ATMOSPHERIC  MODEL  ELEMENTS  ARE  INITIALIZED,  x 

X  FINALLY,  SYSTEM  CONTROL  CONSTANTS  ARE  SET  TO  GENERATE  THE  OUTPUT.  A  X 

X  FURTHER  DESCRIPTION  OF  THE  ABOVE  FUNCTIONS  IS  CONTAINED  IN  THE      x 

X  DESCRIPTION  OF  ARGUMENTS.                                         x 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  X 

X  DESCRIPTION  OF  ARGUMENTS  » 

X  X 

X  INPUT  ARGUMENTS  X 

X  X 

X  PROGRAM  CONTROLS  X 

X  X 

X        ITERl INTEGER,  MAXIMUM  NUMBER  OF  REVOLUTIONS  DESIRED  x 

X  FOR  CONSIDERATION.  A  CONTROL  PARAMETER.  X 

X         ITER2 INTEGER,  ORBITAL  PRINT  FREQUENCY.  A  CONTROL  X 

X  PARAMATER.  SETTING  ITER2  =  TO  100  WILL  REDUCE  X 

X  TABULAR  OUTPUT  BY  ONLY  PRINTING  OUTPUT  EVERY  100  X 

X  ORBITS.  X 

X        PSENS INTEGER,  CONTROL  PARAMETER.  PSENS=1  SPECIFIES  X 

X  SENSITIVITY  ANALYAIS  MODE.  PSENS=0  SPECIFIES  X 

X  STANDARD  PROPELLENT  LONGIVITY  MODE.  X 

X        JET INTEGER,  CONTROL  PARAMETER.  IET=1  PRINTS  THE  X 

X  INITIAL  ANOMALISTIC  PERIOD.  IET=0  SUPRESSES  PRINT. x 

X-  SYSTEM  RESET  TO  0  WHEN  SENSITIVITY  MODE  IS  SELECTED 

X        NOA INTEGER,  NUMBER  OF  SCHEDULED  ORBIT  ADJUSTS.  THIS  X 

X  IS  A  CONTROL  PARAMETER  WITH  ALLOWED  VALUES  FROM  X 

X  ZERO  TO  TEN.   THE  ONLY  ORBIT  ADJUSTS  CURRENTLY  x 

X-  SUPPORTED  ARE  CHANGES  TO  THE  SEMI-MAJOR  AXIS.  X 

X         IRV(K) INTEGER  ARRAY.  (K=10  MAX). THE  ELEMENTS  OF  THE  X 

X  ARRAY  SPECIFY  THE  REVOLUTION  NUMBER  DURING  WHICH  X 

X  EACH  OF  THE  ORBITAL  ADJUSTS  ARE  PERFORMED.  K  MUST  x 

X  EQUAL  NOA.  X 

X        DA(K) INTEGER  ARRAY.  (K=10  MAX)  THE  ELEMENTS  OF  THE  X 

X  ARRAY  SPECIFY  THE  CHANGE  IN  KILOMETERS  OF  THE  X 

X  SEMIMAJOR  AXIS.  K   MUST  EQUAL  NOA.  X 

X        PSIZE INTEGER, CONTROL  PARAMETER  THE  NUMBER  OF  INTERVALS  X 

X  IN  THE  SENSITIVITY  ANALYSIS.  X 

X  X 

X  •    SATELLITE  PHYSICAL  PARAMETERS  X 

X  X 

X        SPI REAL,  MOTOR  SPECIFIC  IMPULSE  IN  SECONDS.  AGRAWAL  X 

X  (REF.  2:PP.  163-166)  GIVES  SAMPLE  VALUES.  X 

X        WT REAL,  INITIAL  MASS  OF  THE  SATELLITE  IN  KILOGRAMS.  X 

X  (INCLUDING  THE  FUEL  MASS)  X 

X        W REAL,  INITIAL  MASS  OF  FUEL  IN  KILOGRAMS.  X 

X         CD REAL,  DRAG  COEFFICIENT  OF  THE  SATELLITE.  X 

X  DIMENSIONLESS.  NORMALLY  VARIES  FROM   ONE  TO  X 

X  AROUND  THREE,  DEPENDING  ON  VELOCITY,  DENSITY,  X 

X  COMPOSITION  OF  ATMOSPHERE  AND  PHYSICAL  DESIGN.  X 

X  IS  MOST  OFTEN  SET  AT  2.20   FOR  ALTITUDES  BETWEEN  X 

X  200  AND  ^00  KM.  FOR  ALTITUDES  OF  1000  KM  CD  IS  x 

X  SET  AT  JUST  ABOVE  3.0.  JACCHIA  (REF.  3:P  2).  x 

X        A REAL,  SATELLITE  CROSSESTIONAL  AREA  IN  SQUARE  X 

X  KILOMETERS.  X 

X  '  X 

X  ATMOSPHERIC  SPECIFIC  PARAMETERS  X 

X  X 

X        D REAL,  RATIO  OF  ATMOSPHERIC  TO  EARTH  ANGULAR  x 

X  ROTATION  RATES.  A  CONSTANT  IN   EQUATION  11  OF  X 

X  NSWC  TR  83-243.  UNITLESS.  THE  ACTUAL  VALUE  WILL  x 

X  DEPEND  ON  LATTITUDE  AND  ALTITUDE.  RANGE  .8  TO  1.0.x 

X         F107 REAL,  DECIMETRIC  SOLAR  FLUX.  JACCHIA  (REF  <+  S  5 )  .  X 

X  THIS  VALUE  IS  USED  BY  THE  INDUSTRY  TO  REPRESENT  x 

X  THE  XUV  FLUX  IN  THE  THERMOSPHERE.  THE  XUV  FLUX  X 
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X  IS  THE  FACTOR  CREATING  THE  SUN  FACING  DIURNAL     x 

»  DENSITY  AND  TEMPERATURE  BULGE  IN  THE  THERMOSPHERE .  5< 

X  IT  IS  DIRECTLY  RELATED  TO  SOLAR  ACTIVITY.   IT      ^ 

*  VARIES  WITH  THE  27  DAY  SLOAR  ROTATION  AND  11  YEAR  3< 

K  SUNSPOT  ACTIVITY  CYCLE.  ORR  (REF.  6:APPENDIX  B)    X 

X  INDICATES  RANGES  OF  73.3  TO  2<42.9.  XUV  FLUX  AND    X 

X  THE  GEOMAGNETIC  HEATING  EFFECTS  CAN  CAUSE          * 

X  DENSITY  TO  CHANGE  BY  THREE  ORDERS  OF  MAGNITUDE  AT  x 

X  600  KILOMETERS.  THE  10.7  CM  FLUX  IS  DECTABLE  FROM  3< 

X  THE  GROUND  AND  THEORETICALLY  PARELLELS  THE  XUV     x 

X  FLUX.  THE  XUV  FLUX  IS  NOT  MEASUREABLE  FROM  THE    x 

X  GROUND  DUE  TO  NEAR  COMPLETE  UPPER  ATMOSPHERE      X 

X  ABSORPTION.                                       X 

X        FBAR REAL,  THE  AVERAGE  10.7  CM  FLUX.  THE  EFFECTS  OF  THEX 

X  XUV  HEATING  DO  NOT  RETURN  TO  THE  ORIGINAL  UNHEATEDX 

X  CONDITION  IN  AN  EARTH  ROTATION.  A  HIGH  HEATING    x 

X  CONDITION  WILL  HAVE  SOME  LONGER  TERM  EFFECTS.      X 

X-  THESE  LONGER  TERM  EFFECTS  ARE  MODELED  HERE  BY     x 

•X  TAKING  A  MOVING  90  DAY  10.7  CM  AVERAGE.            X 

X  FBAR  RANGES  FROM  A  LOW  OF  73  TO  A  HIGH  OF  230.     x 

X        AKP REAL,  GEOMAGNETIC  ACTIVITY  INDEX,  JACCHIA  (REF.  ^    X 

X  AND  5).  GEOMAGNETIC  ACTIVITY  ALSO  CAUSES  HEATING   X 

X  IN  THE  THERMOSPHERE.  THE  EXACT  MECHANISM  IS  NOT    x 

X  UNDERSTOOD.  IDEALLY  QUEIT  MAGNETIC  CONDITIONS      x 

X  CORRESPOND  TO  AN  INDEX  OF  0.  MAXIMUM  ACTIVITY  FOR  X 

X  WHICH  SIGNIFICANT  DATA  EXISTS  ARE  RELATED  TO  AN    X 

X  INDEX  OF  BETWEEN  6  AND  7.  NORMAL  ACTIVITY         » 

X  INDICATES  AN  INDEX  OF  BETWEEN  l.AND  2.             X 

X        TVE REAL,  TIME  OF  VERNAL  EQUINOX  PASSAGE.  USED  TO'     X 

X  TRANSFORM  TIME  FROM  INERTIAL  TO  EARTH  FIXED        X 

X  REFERENCE  FRAMES.  FOUND  IN  THE  ASTRONOMICAL  ALMANAC 

X  FOR  THE  YEAR  OF  CHOICE.  (MODIFIED  JULIAN  DAYS).    x 

X  A  DESCRIPTION  IS  IN  TAFF  (REF.  7:PP.  103-lG<i).     X 

X  X 

X  NAVSPASUR  ORBITAL  ELEMENTS                                  X 

X  X 

X  THE  NAVSPASUR  ORBITAL  ELEMENT  SET  IS  THE  BROUWER  MEAN       X 

X  ELEMENT  SET.  BROUWER  (REF.  1)  IS  THE  BEST  DESCRIPTION  FOR    x 

X  THIS  ELEMENT  SET.                                         X 

X  X 

X        UJD REAL,  TIME  OF  EPOCH  IN  MODIFIED  JULIAN  DAYS.       X 

X  THIS  IS  THE  TO  IN  THE  EQUATION  FOR  THE  MEAN  ANOMALY 

X  M.  (  M  =  MO  +  N  X  (T  -  TO)  ) .  TAFF  (REF.  7).       X 

X         ETU REAL,  INITIAL  MEAN  ANOMALY  MO.  MEASURED  IN  DEGREESX 

X  HAS  A  RANGE  OF  0  TO  360.                          X 

X        GO REAL,  MEAN  ARGUMENT  OF  PERIGEE.  RANGES  FROM  0  TO   x 

X  360  DEGREES.  NOT  DEFINED  FOR  CIRCULAR  ORBITS.      x 

X        HO REAL,  MEAN  RIGHT  ASCENSION  OF  THE  ASCENDING  NODE.  X 

X  RANGES  FROM  0  TO  360  DEGREES.  NOT  DEFINED  FOR     x 

X  EQUITORIAL  ORBITS                                 X 

X        ES REAL,  MEAN  ECCENTRIY.  FOR  THIS  THEORY  RANGES  FROM  x 

X  0.01  TO  0.1.  THESE  ARE  LOW  ECCENTRICITY  ORBITS.    X 

X        XI REAL,  MEAN  INCLINATION.  RANGES  FROM  0  TO  180       X 

X  DEGREES.  VALUES  OF  0  AND  130  ARE  IN  EQUITORIAL     x 

X  ORBITS.  INCLINATION  OF  90  IS  A  POLAR  ORBIT.  RANGESX 

X  -        OF  0  TO  90  ARE  IN  PROGRADE  ORBITS.  INCLINATIONS    X 

X  90  TO  180  ARE  IN  RETROGRADE  ORBITS.  FOR  Afl         X 

X  INCLINATION  OF  0  OR  180  HO  IS  NOT  DEFINED.  THERE   X 

X  IS  NO  ASCENDING  NODE  IN  EQUITORIAL  ORBITS.         X 

X        RND REAL,  RATE  OF  CHANGE  OF  MEAN  MOTION.  MEASURED  IN   X 

X  REVOLUTIONS  PER  (DAY  SQUARED) . VALUES  OF  .00001     X 

X  AND  ABOVE  ARE  CONSIDERED  SIGNIFICANT  WHEN          X 

X  ATMOSPHERE  INDUCED.  REFERENCE  (TAFF)      P137      X 

X        ESD REAL,  ECCENTRICITY  DECAY  RATE  MEASURED  IN         X 

X  INVERSE  SECONDS.                                 x 

X        AO REAL,  KAULA  SEMIMAJOR  AXIS.  MEASURED  IN  X 

X  MEAN  EQUITORIAL  EARTH  RADII  (6378.165  KM).  RANGES  x 

X  FROM  1.02594  TO  1.168^+7  TO  MAINTAIN  AN  ALTITUDE  OFX 

X  100  TO  1000  KM  FOR  RADIUS  AT  PERIGEE.  THIS  MUST  BEX 

X  CAREFULLY  SCALED  TO  ECCENTRICITY  TO  PREVENT        X 

X  CRASHING  THE  SATELLITE.  PERIGEE  ALTITUDES  ABOVE   x 
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X                  600  KM  ARE  ESCENTIALLY  AVOVE  THE  ATMOSPHERE.  X 

X        ADOT REAL,  SEMIMAJOR  AXIS  DECAY  RATE.  X 

X  X 

X  X 

X  K 

X  X 
XXXXXXXX?f3<XXXX3<XXXXXXXXXXXXXXX3(3(50(K3(X-«XXXXXXXXXXXXXXXXX»9«XKXXXXXXXXKXX 
X 

DIMENSION  SNSAR(20),  SNSMIN(2Q),  IRV(IO),  DA(IO),  SN5DEL(20) 
X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X    THE  SET  OADJ  IS  COMMON  TO  PLEP,  KOZAI,  THRST  AND  SNALYS.  X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

COMMON  /  OADJ/  NOA,  IRV,  DA,  SPI,  WT,  PSENS 
X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X    THE  SET  XXXX  IS  COMMON  TO  PLEP  AND  PERIOD.  x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

COMMON  /  XXXX  /  lET 
X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X    THE  SET  SOLAR  IS  COMMON  TO  PLEP  AND  DRAG.  X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

COMMON  /  SOLAR  /  F107  ,  FBAR  ,  AlCP  ,  TVE 
INTEGER  PSIZE,  PSENS,  PSYNS 
REAL  MINSN 
X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X    INITIALIZING  CONSTANTS.  X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

DATA    R2   /    541.15E-06    /    ,DEG  /57 .295779513D0   /    ,    ROE  /6375.165D0/ 
X 

xxxxxxx=<xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X         INPUTTING   CONTROL    PARAMETERS.  X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

READ  X,  ITERl,  ITER2,  NOA,  PSENS,  lET,  PSIZE 

X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X    INPUTTING  SATELITE  SPECIFIC  PARAMETERS.  X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

READ  X,  CD,  A,  D,  SPI,  W,  WT 

X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X    INPUTTING  EARTH  AND  ATMOSPHERIC  PARAMETERS.  x 

XXXX'xXXXXXXXXXXJtXXKXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
X 

READ  X  ,  F107  ,  FBAR  ,  AKP  ,  TVE 

X 

XXXXXX=<^X5(^XXXXX?(XX-^XXXXXXXXXXXX5(XXXXXXXXXX5<^XXXXXXXXXXXXXXXXXJ(XX?(XXXX 
X         INITIALIZING    BROUWER    MEAN    ORBITAL    ELEMENTS    FROM    THE    NAVSPASUR  ^ 

X         FIVE    CARD    FORMAT    ORBITAL    DATA    SET.  X 

KXXXxXXXXXXX?<XXXXXXXXX;;XXXXXX5?XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
X 

READ    5    ,    UJD    ,    ETU    ,    HO    ,    GO    ,    ES    ,    XI 
5    FORMAT    (    /8X,F14.8,5(lX,F8.<i)) 
READ    10    ,    RND    ,    ESD 
10    FORMAT    (    20X,    F11.9,    19X,E12.5    ) 

READ    15    ,    AO    ,    ADOT 
15    FORMAT    (    /8X,F11.5,1X,E14.7) 
X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X         INPUTTING   THE   REVOLUTION    NUMBER   AND   CORRESPONDING   CHANGE    IN  x 

X        SEMIMAJOR    AXIS  X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

IF    (    NOA    .EQ.    0)    GO    TO    30 
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DO    20      K    =    1,    NOA 

READ    X,    IRV(K)    ,    DA(K) 
20    CONTINUE 

XXXXXXXXX5<5«XXXX3(XX5(XX5(9<XXXXX=<XXXX3(3(XXX3(5(XX5<X5<XXX5<XX3<X3(XXXX5(X5«X5(3<XKXXXX 

X  THIS  NEXT  SERIES  OF  CODE  SETS  THE  INITIAL  VALUES  FOR  THE  SCANNING  X 
X  ARRAY  MATRIX  SO  AS  TO  ALLOW  THE  RECOVERY  OF  THE  ORIGINAL  INPUT  X 
X   ELEMENTS  IF  SENSITIVITY  MODE  IS  NOT  SELECTED.  x 

XXXXXXXXXXXX5(XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX3(XXXXXXXXXXXXXXXXXXXXXXX 

30  DO  90  ICN  =  1,20 

SNSDEL(ICN)  =  1.0 

SNSMIN(ICN)  =  0.0 

SNSAR(ICN)  =  0.0 
90  CONTINUE 

XXXXXXXXXXX5(XX5<XXXXXXXXXXXXXXXX3(XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX5(XXXXXX 

X  THE  NEXT  SERIES  OF  WRITE  STATEMENTS  OUTPUTS  THE  ORIGINAL  INPUT      X 

•X  PARAMETERS.  UNITS  ARE  SUPPLIED  FOR  CLARIFICATION.  X 

XXXXX3<XXXXXXXXXXXX3(XXXXX3(X3(XX3<XXXXXXXXXXXXKXXXX5(XXXXXXXXXXXXXXXXXXXXXX 

WRITE(6,98) 
WRITE(6,99) 
WRITE(6,103) 
WRITEC6,99) 

WRITE(6,10<4)  D,  F107,  FBAR,  AKP,  TVE 
WRITE(6,105)  CD,  A,  SPI,  W,  WT 
WRITE(6,106)  UJD,  ETU,  HO,  GO,  XI 
WRITE(6,107)  RND,  ES,  ESD,  AO,  ADOT 
WRITE(6,10S)  ITERl,  NOA 
IFCNOA  .GT.  0)  THEN 
WRITE(6,109) 
DO  27  NN  =  1,  NOA 

WRITE(6,110)  IRV(NN),  DA(NN) 
27       CONTINUE 
END  IF 
WRITE(6,99) 

XXXXXXXX5(XXXXXXXXXX9(XXXX?(X3(X3«5(XXXXXXXXXXXXXXXXXXXXXXX5(XXXX5(XXX5(5<XXXXXX 

X  THIS  NEXT  CODE  SETS  THE  INTIAL  INPUT  PARAMETERS  TO  SOME  RETAINERX 
X  VARIABLES  THAT  ALLOW  RETURNING  TO  THE  ORIGINAL  VALUES  IF  THE  X 
X  SENSITIVITY  MODE  IS  SELECTED.  THIS  PREVENTS  HAVING  TO  REREAD  ALL  X 
X  THE  INPUT  PARAMETERS  BACK  IN  FOR  EACH  RUN  OF  THE  SENSITIVITY  SCAN.X 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

V1=CD 

V2=A 

V3  =  D 

V4=SPI 

V5=W 

V6=WT 

V7=F107 

V3=FBAR 

V9=AKP 

V10=TVE 

V11=UJD 

V12=ETU 

V15=H0 

Vl'i  =  GO 

V15=ES 

V16=XI 

V17=RND 

V13=ESD 

V19=A0 

V20=ADOT 

XX3(XXX3(XXXXX30(XXXXX5<X5<XX5(XXXXXXXX3(XXX5(XXX)(XXXXXX3(XXXX50<XXXXXX5(XXXXXXXX 
X  THIS  IS  THE  ITERATIVE  SCANNING  ROUTINE.  ARGUMENTS  ARE  INCREMENTED  X 
X  FRACTIONALLY  THROUGH  THEIR  SPECIFIED  RANGE  AS  SPECIFIED  IN  THE  X 
X  INPUT  TO  THE  SUBROUTINE  SNALYS  (P5IZE  ITERATIONS  OF  PLEP).  X 
X   IF  SENSITIVITY  MODE  IS  NOT  SELECTED  PSIZE=1.  X 

XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXJ(XXXXXXXXXXXXXXXXXXXXXXXXXXX3(XXXXXX 
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IFCPSENS  .EQ.  0)  THEN 

PSIZE=0 
END  IF 

DO  ^0  PSYNS  =  0,  PSI2E 
X 

X  THIS  NEXT  SET  OF  CODE  IS  ENVOKED  WHEN  PSENS  IS  EQUAL  TO  ONE.  x 
X  IT  SETS  UP  FOR  AND  CALLS  THE  SENSITIVITY  ARGUMENT  ARRAY  SNALYS  FORX 
*      SUBSEQUENT  SCANNING.  (PSENS  =  1  SELECTS  SENSITIVITY  MODE  x 

X3(X3(3<K3(XX3(XXXXXXXX5<XXXXXXXK3<XXXX3(X5(XXXXXXXX3()0(XXXX3«X5«3<XXXXXXXXXXXX)(5(XX 

IF  (PSENS  .EQ.  1)  THEN 
IET  =  0 

XX)(XXXKXXXXXX3(XX3(XX5(XXX3(XX5<XXXX3(5(30«XXXXXXXXXXX3(9<XXXXX?(XXXXXXXXXXXXXXXX 

X  THE  SUBROUTINE  SNALYS  INTIALIZES  THE  SCANNED  INPUT  PARAMETERS  TO  x 
X  BE  INCLUDED  IN  THE  STUDY.  IT  ALSO  PRINTS  THE  SENSITIVITY  MATRIX.  x 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

IF(PSYNS  .EQ.  0)  THEN 

CALL  SNALYS(  PSIZE,  SNSAR,  SNSMIN,  SNSDEL,  SUMCK,  DELSN, 
1  MINSN) 

ELSE 

CD  =   VI 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X  RECOVERY  OF  THE  ORIGINAL  INPUT  PARAMETERS  AFTER  THE  FIRST  RUN  OF  x 
X  A  SENSITIVITY  ANALYSIS.  x 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

A   =   V2 

D   =   V3 

SPI=  v^ 

W   =   V5 

WT  =   V6 

F107=  V7 

FBAR=  V8 

AKP=   V9 

TVE=V10 

UJD=V11 

ETU=VI2 

HO  =V13 

GO  =V14 

ES  =V1S 

XI  =V16 

RND=V17 

ESD=V18 

AO  =V19 

ADOT=V20 
END  IF 
END  IF 

ASYNS=REAL(PSYNS) 
ASIZE=REAL(PSIZE) 
IF(PSENS  .EQ.  1)  THEN 

Ir(SUMCK  .EQ.  1.0)  THEN 

SCNFRC  =  ASYNSXDELSN  +  MINSN 
ELSE 

SCNFRC  =  ASYNS/ASIZE 
END  IF 
END  IF 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X  THIS  NEXT  CODE  CREATES  THE  CURRENT  VALUE  FOR  THE  SCANNED  INPUTS.  X 
X  IF  SENSITIVITY  IS  SELECTED  THE  PARAMETER  WILL  SCAN  FROM  MIN  TO  MAXX 
X   WITH  A  GRANULARITY  OF  1/PSIZE  PER  SCAN  ^ 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

CD  =  SNSAR(1)X(SNSMIN(1)+ASYNSXSNSDEL(1))-CDX(SNSAR(1)-I.0) 

A   =  SNSAR(2)X(SNSMIN(2)+ASYNSXSNSDEL(2))-A  x( SNSAR( 2)-l . 0 ) 

D   =  SNSAR(3)X(SNSMIN(3)+ASYNSXSNSDEL(3))-D  x(SNSAR( 3)-l . 0 ) 

SPI=  SNSAR(4)X(SNSMIN(^)+ASYNSXSNSDEL(4))-SPIX(SNSAR('4)-1.0) 
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W   =   SMSAR(5)X(SNSMIN(5)+ASYNS5(SNSDEL(5))-N  X(  SNSAR(  5)-l  .  0  ) 
WT  =   SNSAR(6)X(SNSMIN(6)+ASYNS)(SNSDEL(6))-WTX(SNSAR(6)-1.0) 
F107=  SNSAR(7)x(SNSMim7)+ASYNSXSNSDEL(7))-F1073<(SNSAR(7)-1.0) 
FBAR=  SNSAR(8)X(SNSMIN(8)  +  ASYNSXSNSDEL(3))-FBAR5((SNSAR(8)-1.0) 
AKP=   SNSAR(9)i<(SNSMIN(9)+ASYNSXSNSDEL(9))-AKP^(SMSAR(9)-1.0) 
TVE=SNSAR(10)X(SNSMIN(10)+ASYNS3(SNSDEL(10))-TVE5((SNSAR(10)-1.0) 
UJD  =  SNSARC11))((SNSMIN(11)+ASYNS5(SNSDEL(11))-UJDJ((SNSAR(11)-1.0) 
ETU  =  SNSAR(12  3X(SNSMIN(12)+ASYNSXSNSDEL(12))-ETU5((SNSAR(12)-1.0) 
HO  =SNSAR(13)X(SMSMIN(13)+ASYNSXSNSD£L(13))-H0X(SMSAR(13)-1.0) 
GO  =SNSAR(l<t)X(SMSMINC14)+ASYNS5<SNSDEL(14))-G03((SNSAR(l'^)-l  .0) 
ES  =SNSAR(15)X(SNSMIN(15)+ASYNS^^SNSDEL(15))-ESX(SNSAR(15)-1.0) 
XI  =SNSAR(16)X(SNSMIN(16)+ASYNS5<SMSDEL(16))-XIX(SNSAR(16)-1.0) 
RND  =  SNSAR(17)X(SNSMirU17)+ASYNS5(SNSDEL(17))-RND5((SNSARa7)-1.0) 
ESD  =  SNSAR(18)X(SNSr'1IN(18)  +  ASYNS3(SNSDEL(lS))-ESDX(SNSAR(18)-1.0) 
AO  =SNSAR(19)X(SNSMIH(19)+ASYNSXSNSDEL(19))-A05«(SNSAR(19)-1.0) 
ADOT  =  SMSAR(20)X(SNSMIN(20)+ASYNS5<SNSDEL(20))-ADOTX(SNSAR(20)-1.0) 

X  DEL  IS  THE  CONSTANT  MULTIPLIER  FOR  THE  INTEGRAL  IN  EQUATION  33  OF  X 
X  NSWC  TR  83-2<+3  X 

xxxxxxxxxxx;<xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

DEL  =  CDXAX3.14159D0/SPI  . 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X  IF  THE  SENSITIVITY  ANALYSIS  MODE  IS  NOT  SELECTED,  THEN  THE  SYSTEM  X 
X   PRINTS  STANDARD  MISSION  LIFE  TIME  TABULAR  DATA.  X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

IF  (PSENS  .EQ.  0)  THEN 

BCOF  =  CD  X  A 

WRITE  (6,98) 

WRITE  (6,99) 

WRITE  (6,100) 

WRITE  (6,99) 

WRITE  (6,101)  BCOF,  SPI,  W 

WRITE  (6,102)  F107,  FBAR,  TVE 
END  IF 
X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
X   CONVERTING  DEGREES  TO  RADIANS  X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

XI  =  XI  /  DEG 
GO  =  GO  /  DEG 
HO  =  HO  /  DEG 
FLff  =  ETU  /  DEG 

X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  THIS  NEXT  CODE  SEQUENCE  CONDITIONS  THE  ELEMENTS  AS  PER  NSWC  TR  82-387 
X  PAGE  7.  IT  CONVERTS  THE  KAULA  SEMIMAJOR  AXIS  AO  IN  EARTH  RADII  TO  X 
X  BROUWER  MEAN  SEMIMAJOR  AXIS  VIA  EQNS  2.1  AND  2.2.  THE  SEMIMAJOR  AXISX 
X  DECAY  RATE  ADOT  IS  CONVERTED  IN  A  LIKE  MANNER  VIA  2.3,  2.(4,  AND  2.5  X 
XXXX^XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
X 

FN  =  1.  /  (  1.  -  ES  X  ES  )XX1.5D0 
TA  =  SIN  (  XI  ) 
TB  =  TA  X  TA 
TE  =  1.  -  (3.XTB)  /  2. 
TAD  =  ((  3.XR2)  /  (2.XA0XA0))XTEXFN 
AK  =  AO 

AO  =  (A0X((1.  +  2.XTAD)  /  (1.  -  TAD) )XX0 . 66666667D0 )XROE 
XDOT  =  TADX(3.XESDX(ES/(1.-ESXES))  -  2  .  X( ADOT/AK) ) 
ADT  =  (ADOT/AK)XA0  +  2 . XAKXR0EXXD0TX( ( 1 .  +  2.xTAD)x 
1  ((1.  -  TAD)xx  5))xx(-0.333-3333333D0) 

XXXXXXXXXXXXXXXXXXXXXXXXXXX3<XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

X  TIME  IS  NOW  CONVERTED  FROM  SECONDS  TO  DAYS  X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

ADT  =  ADT  /  86400. DO 
ESD  =  ESD  /  86400. DO 
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RND  =  RND5«2.D0X3.14159D0X(1./(86400.D0)XX2) 

X 
XXX3(XXX50fXX5(X3(XXX3(X^3(XX5(X5(XX5(X5«XXXXX=<XXXXXX5<XXXXXXXXXX)(X)(XXXX3(3(XXXXXX3(X 

X  THE  SUBROUTINE  KOZAI  IS  THE  OPERATING  WORKHORSE  FOR  THE  PROPELLENT  x 
X  LONGEVITY  COMPUTATIONS.  IT  ALSO  DECREMENTS  PROPELLENT  AS  IT  IS  USEDX 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

CALL  KOZAI(ADT,AO,ES,XI,FLO,GO,HO,UJD,W,DEL,ITER1,ITER2, 
1  SCNFRC,  D) 

X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X   FORMAT  STATEMENTS  USED  IN  THE  PLEP  INITIALIZATION  PROGRAM.         X 

XXXXXXXXXXXX?{XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
X 

98  FORMAT  (  IHl  ) 

99  F0RMAT(//5X,  ' ',//) 

100  F0RMAT(//I5X, 'PROPELLANT  LIFE  ESTIMATOR',//) 

101  F0RMAT(//5X, 'CDA=',E12.5, •  SQ  KM'/5X, 'SPECIFIC  IMPULSE=', 
XF7.2,'  SEC'/5X, 'INITIAL  WEIGHT  OF  PROPELLANT= • , 

XF7.2,'  KG') 

102  FORMAT (//5X, ' F107= ' , F8 . 2/5X, ' FBAR= ' , F8 . 2/5X, 
X  '  TIME  OF  VERNAL  EQUINOX=  ' ,  E22  .  l-i  ) 

103  FORMATdSX,  'INITIAL  INPUT  CONDITIONS') 

104  F0RMAT(///1X, 'ATMOSPHERIC  INPUT  PARAMETERS '//5X, ' D  (RELATIVE  ATMOS 
XPHERIC  ROTATION)  =  ' ,  FIO .  <4/5X,  '  F107  (DECIMETRIC  SOLAR  FLUX)  =  ', 
XF7.2/5X, 'FBAR  (90  DAY  AVG.  F107)  =  ' , F7 . 2/5X, ' AKP  (GEOMAGNETIC  IND 
XEX)  =  •,F6.3/5X, 'TVE  (TIME  OF  VERNAL  EQUINOX  PASSAGE)  =  ',F16.9, 
X'  MODIFIED  JULIAN  DAYS'///) 

105  FORMATdX, 'SATELLITE  PHYSICAL  PARAMETERS  «//5X, 'CD  (DRAG  COEFFICIEN 
XT)  =  ',F5.2,'  UNITLESS'/5X, 'A  (CROSSECTIONAL  AREA)  =  ',F11.9,'  SQU 
XARE  KIL0METERS'/5X, 'SPI  (MOTOR  SPECIFIC  IMPULSE)  =  ',F8.2,'  SECOND 
XS'/5X,'W  (INITIAL  FUEL  MASS)  =  •,F9.2,'  KIL0GRAMS'/5X, 'WT  (INITIAL 
X  SATELLITE  MASS)  =  '.F11.2,'  KILOGRAMS'///) 

106  FORMAT(  IX, ' NAVSPASUR  ORBITAL  ELEMENTS '//5X, ' UJD  (EPOCH)  =  •,F14.8 
X,'  MEAN  JULIAN  DAYS '/5X, ' ETU  (MEAN  ANOMALY)  =  '^FS.A,'  DEGREES'/5X 
X,'HO  (MEAN  RIGHT  ASCENSION  OF  THE  ASCENDING  NODE)  =  ',F8.4,'  DEGRE 
XES'/5X,'G0  (MEAN  ARGUMENT  OF  PERIGEE)  =  ',F8,4,'  DEGREES'/5X, 
X'XI  (MEAN  INCLINATION)  =  •,F8.<4,'  DEGREES') 

107  F0RMAT(5X, 'RND  (RATE  OF  CHANGE  OF  MEAN  MOTION)  =  SFll.g,'  REVOLUT 
XTIONS  PER  DAY  SQUARED'/ 

X5X,'ES  (MEAN  ECCENTRICITY)  =  •,F8.4,'  UNITLESS'/5X, • ESD  (ECCENTRIC 
XITY  DECAY  RATE)  =  •,F16.9,*  1/SECONDS '/5X, ' AO  (KAULA  SEMI-MAJOR  AX 
XIS)  =  ',F11.5,'  MEAN  EAKTH  RADII '/SX, 'ADOT  (SEMI-MAJOR  AXIS  DECAY 
XRATE)  =  ',F16.9,'         '///) 

108  FORMATdX, 'CONTROL  INPUT  PARAMETERS'//5X, 'MAXIMUM  REVOLUTIONS  PER 

.  XITERATION  ', 17 ,/5X, ' NUMBER  OF  SCHEDULED  ORBIT  ADJUSTS  TO  SEMI-MAJO 
XR  AXIS  '  12/ ) 

109  F0RMAT(//5X, 'SEMI-MAJOR  AXIS  INCREASE  SCHEDULE'//5X, 'REV  ADJUST', 
X5X, 'CHANGE  IN  KILOMETERS'/) 

110  F0RMAT(5X,I7,7X,F6.2) 


40  CONTINUE 

::ap 

END 
X 
X 
X 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  XX  XX 

X  XX    SUBROUTINE  SNALYS     xx 

X  XX  .  XX 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 
X 

X 

SUBROUTINE  SNALYS  (PSIZE, SNSAR,  SNSMIN,  SNSDEL,  SUMCK,  DELSN, 
IMINSN) 
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DIMENSION  SNSAR(20),  SHSMIN(20),  SNSMAX(20),  SNSDEL(20) 

CHARACTER5<A  SNSCAR 

INTEGER  N,  PSIZE 

REAL  SNSAR,  SNSMIN,  SNSMAX,  SHSDEL,  SUMCK,  DELSN,  MINSN 

WRITE  (6,98) 

WRITE(6,120) 

WRITE(6,100) 

WRITE(6,120) 

WRITE(6,200)  PSIZE 

ASIZE=REAL(PSIZE) 

READ(5,102) 

SUMCX  =0.0 

DO  20  N  =  1,  20 

READ  101,  SNSCAR,  SNSAR(N),  SNSMIN(N),  SNSMAX(N) 
SNSDEL(N)  =  (SNSMAX(N)-SNSMIN(N))/ASIZE 
IF  (SNSAR(N)  .EQ.  1.0)  THEN 
SUMCX  =  SUMCK  +1.0 
DELSN  =  SNSDEL(N) 
MINSN  =  SNSMIN(N) 

WRITE(6,300)  SNSCAR,  SNSMIN(N),  SNSMAX(N), 
1  SNSDEL(N) 

END  IF 
20  CONTINUE 

IF  (SUMCK  .GT.  1.0)  THEN 
DELSN  =  0.0 
MINSN  =0.0 
END  IF 

WRITE(6,120) 
WRITE(6,^00) 

98  FORMAT  (  IHl  ) 
120  F0RMAT(//5X,  ' • ',//) 

100  F0RMAT(//17X, 'SENSITIVITY  ANALYSIS  FOR ' ,/16X, ' LOW  EARTH  ORBIT  ', 
X  'SATELLITES', /19X, 'PROPELLENT  LONGEVITY',//) 

200  FORMAT  (5X,  '  TABLE  OF  INPUT  PARAMETERS  SCANNED  FROM  MINIMUM  TO', 
X  /5X,'  MAXIMUM  WITH  THE  INDICATED  INCREMENT  PER  RUN.  ',/5X, 
X  •  THE  GRANULARITY  IS  1/ ' , 14, ' . '/5X,  •  ALL  OTHER  INPUT  PARAMETERS 
X  ARE  RETURNED  TO  THEIR', /5X,'  INITIAL  VALUES  AT  EACH  INCREMENT.' 
X  ///'INPUT  PARAMETER' ,8X, 'MINIMUM', 13X, 'MAXIMUM' , 13X, 'INCREMENT'/) 

102  FORMAT  (IX) 

101  FORMAT  (  A-^,  IX,  F3.1,  F16.9,  4X,  F16.9) 
300  FORMAT  (5X,  A4,  6X,  3(3X,  F16.9)) 

<»00  FORMAT(//,  lOX,  '  SENSITIVITY  SCAN  RESULTS  FOR  ABOVE  SPECIFICS'// 
15X,'REV  NUMBER', 8X, 'TIME(DAYS)',15X, 'FRACTION  OF  SCAN') 

RETURN 
END 


X 
X 

X 
X 

X  XXXXXKXXXXXX»XX5<XX3<XXXXXXXXXXX 

X  XX  XX 

X  XX      SUBROUTINE  KOZAI      xx 

X  XX  XX 
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X  XXXX5(X5<J<5(XX3«XXXXXX3(XXX5(XXXXXXX 

X 

X 

X  X 

X   SUBROUTINE  KOZAI  IS  THE  SYSTEM  WORKHOUSE.  IT  CALLS  BRAUER  TO  CON-x 

X  VERT  THE  BROUWER  MEAN  ORBITAL  ELEMENTS  OF  THE  NAVSPASUR  INPUT  TO  x 

X  BROUWER  OSCULATING  ELEMENTS.   THESE  OSCULATING  ELEMENTS  ARE  CON-  x 

X  VERTED  TO  KOZAI-LIKE  MEAN  ELEMENTS  IN  THE  SUBROUTINE  MEAN.  THE  x 

X  SUBROUTINE  PERIOD  COMPUTES  THE  INITIAL  ANOMALISTIC  PERIOD.  THE  x 

X  SUBROUTINE  GEOP  COMPUTES  THE  EFFECTS  OF  THE  OBLATE  GEOID  ON  THE  X 

X  ORBITAL  ELEMENTS  (THROUGH  THE  J4  TERM).  SUBROUTINE  THRST  COMPUTES  x 

X  THE  EFFECTS  ON  THE  ELEMENTS  DUE  TO   PERIGEE  BURNS  TO  CHANGE  X 

X  THE  SEMIMAJOR  AXIS  DURING  IRV(K)  BY  THE  AMOUNT  DA(K)  .  FINALLY,  x 

X  SUBROUTINE  DRAG  COMPUTES  THE  THE  AMOUNT  OF  FUEL  NEEDED  TO  OVERCOME  X 

X  THE  DRAG  EFFECTS  BY  CONTINUOUS  INTRACK- MICROTHRUSTING.  x 

X  KOZAI  KEEPS  TRACK  OF  THE  FUEL  USED   AND  ORBITAL  ELEMENT  CHANGES.  X 

X  WHEN  THE  FUEL  IS  EXHAUSTED  KOZAI  RETURNS  SYSTEM  CONTROL  TO  PLEP.  X 

X  X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 
X 

SUBROUTINE  KOZAKADT,  AO  ,  ES,  XI,  FLO,GO,  HO,  UJD,  W,  DEL,  ITERl,  ITER2, 
1  SCNFRC,  D) 

DIMENSION  X0SC(6)  ,    XM(6),  DXDT(6)  ,  XDD(6),  IRV(IO),  DA(IO) 

DIMENSION  DXM(6) 

COMMON  /  OADJ/  NOA,  IRV,  DA,  SPI,  WT,  PSENS 

DATA  B0,B2,B3,B4,B5,B6  /  . SgseOSZS^E+Ofi , - . 1755528999E+11, 
1      .26386647738E+12, . 1063073996E+16, .805605022E+18, 
1      .7292I15856E-0<i/ 

DATA  DEG/57. 29577951 3D0/,CN2/0.0/,Al, El, RN1/0.0,0.0,0.0/,DT/0.0/ 

DATA  PI  /  S.l^lSgDO  / 

DATA  NPLOT  /  0  / 

IRC  =  0  - 

PCOUNT=0 

KK  =  0 

PI2  =  2.D0X  PI 

CALL  BRAUER  ( BO, B2, B3, B4, B5, DT, AO, ES,XI, FLO,GO, H0,CN2, A,CE,CI,CL , 
1     G,H,A1,E1,RN1  ) 

IF  (PSENS. EQ.  0)  THEN 

WRITE(6,100)  ITERl 

100  F0RMAT(//5X, 'MAXIMUM  NUMBER  OF  ORBITAL  REVOLUTI0NS= • , 110) 
WRITE(6,101)  ITER2 

101  F0RMAT(//5X, 'DATA  PRINTED  EVERY  M3,  •  ORBITAL  REVOLUTIONS') 
SCNFRC  =0.0 

END  IF 

CALL  PERIOD  (  AO,  ES,  XI,  PA  ) 

XOSC(l)  =  A 

X0SC(2)  =  CE 

X0SC(3)  =  CI 

XQSC(^)  =  G 

X05C(5)  =  H 

X0SC(6)  =  CL 

CALL  MEAN  (XOSCXM) 

ILINE  =0 

CIW  =  XM(3)XDEG 

GW  =  XM(4)XDEG 

HW  =  XM(5)XDEG 

CLW  =  XM(6)XDEG 

IF  (PSENS  .EQ.  0)  THEN 
WRITE(6,102) 
WRITE(6,103)  (XM(K),K=1,2),CIW,GW,HW,CLW 

END  IF 
103  F0RMAT(/5X, 'SEMIMAJOR  AXIS=  • ,  E22  .  1<+,  '  KM'/ 
X  5X, •ECCENTRICITY=',E22.14/ 
X  5X, 'INCLINATICN=',E22.1<i,  '  DEG'/ 
X  5X, 'ARGUMENT  OF  PERIGEE= ' , E22 . 14, '  DEG'/ 
X  5X, 'ASCENDING  NODE=  '  ,  E22  .  l<i,  '  DEG'/ 
X  5X,'MEAN  AN0MALY=',E22.14, '  DEG') 
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102  F0RMAT(//5X, 'KOZAI  MEAN  ELEMENT  SET') 
DELT  =  PA 

CALL  THRST  (  XM,  DXM,  IRC,  KK,  PA,  DWOA) 
DO  10    1=1,  ITERl 

T  =  (  I  -  1  )  X  DELT  /  86400. 

TT  =  UJD  +  T 

15  CALL  GEOP  (XM,DXDT) 

16  DO  20    JJ  =  1,6 

XM(JJ)  =  XM(JJ)  +  DXDT(JJ)XDELT  +  DXM(JJ) 
IF  (  JJ  .  LT  .  <+  )  GO  TO  20 

XM(  JJ  )  =  AMOD  (  XM(JJ),  PI2  ) 
20       CONTINUE 
IRC  =  0 
DO  30   II  =  1,  NOA 

IF  (  I.EQ.  IRV(II))  IRC  =  1 
IF  (  IRC  .EQ.  1)  KK  =  II 
30      CONTINUE 

CALL  THRST(XM,  DXM,  IRC,  KK,  PA,  DWOA) 
K 

X    CALCULATING    THE    DELI    IN    NSWC   TR    83-243    EQN    11.  X 

XX5<XXXXXXX5(XXX=<:<XXXXX3(5(XXXXX3(XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX5(XXXXXXX 
X 

RP=XM(1)X(1.0-XM(2)) 

VP=SQRT(CB0/XM(1))X((1.0+XM(2))/(1.0-XM(2)))) 

DEL1=(1.0-(RP/VP)XDXB6XCOS(XM(3)))»X2 

DDELT=DELXDEL1 

CALL  DRAG  (  XM,  DDELT,  TT,  DELMG  ) 
W  =  W  +  DELMG  -  DWOA 
NT  =  WT  +  DELMG  -  DWOA 
IF  (PSENS  .EQ.  0)  THEN 

IF  (PC0UNT.NE.(ITER2-1))  THEN 
PC0UNT=PC0UNT+1 
GO  TO  10 
ELSE 

PCOUNT=0 
END  IF 

IF  (  ILINE  .GE.  50  )  ILINE  =  0 
ILINE  =  ILINE  +  1 

IF  (  ILINE  .EQ.  1  )  WRITE  (  6,104) 
IF  (  ILINE  .EQ.  1)  WRITE  (  6,105) 
WRITE  (  6,  106)  I,  TT,  T,  W 
END  IF 
IF  (  W.LE.  0.00)  THEN 

IF  (NPLOT  .EQ.  1)  THEN 

PRINT  X,  T,  SCNFRC 
ELSE 

WRITE(6,107)  I,  T,  SCNFRC 
END  IF 
RETURN 
END  IF 
10  CONTINUE 

IF  (NPLOT  .EQ.  1)  THEN 

PRINT  X,  T,  SCNFRC 
ELSE 

WRITE(6,107)  I, T, SCNFRC 
END  IF 

104  FORMAT(lHl) 

105  F0RMAT(/5X, 'REV  NUMBER ', 3X, ' TIME(MJD) ', 8X, 'TIMEC DAYS) ', 8X, 
X  'PROPELLANT  REMAININGC KG) ' ,// ) 

106  FORMAT(5X,I10,3X,F9.2,3X,F10.2,15X,F8.2) 

107  FORMAT(5X,I10,  8X, FIO . 2, 15X, F16 . 9 ) 
RETURN 

END 
X 
X 

X  xxxxxxxxxxxxxxxxxxxxx-«xxxxxxxx 

X  xxxxxxxx^xxxxxxxxxxxxxxxxxxxxx 

X  XX  XX 

X  XX    SUBROUTINE  BRAUER     XX 

X  XX  XX 
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X  XXXXXXX5<X)(XX=<XJ(X5(3f5<X3()(y;XXXXX3(K 

K 

X 

SUBROUTINE    BRAUERCBO, B2, B3, B4, B5, DT, A2P, E2P,CI2P, CL02P,G02P, H02P, 
1CN2,A,CE,CI.CL,G,H,ADT,RND,ESD) 
1004   A2P2  =  A2P5«M2 

A2P4  =  A2P2X5(2 

CN0  =  SQRT(B0/(A2P25<A2P)) 

c2P2  =  EZP5«5(2 

ETA=SQRT(1.-E2P2) 

SINEI=SIN(CI2P) 

THETA=C0S(CI2P) 

THETA2  =  THETA5(3(2 

THETA^  =  THETA2^5(2 

THETA6=THETA4*THETA2 

CJ2  =  -B2/(2.?(B0XA2P2) 

ETA2=ETAXX2 

ETA3  =  ETA25<ETA 

ETA<+  =  ETA25t5(2 

CJ21P  =  CJ2/ETA<4 

CJ31P  =  33/(B0XA2P25«A2P5(ETA4)(ETA2) 

CJ'ilP  =  (3.XB<;)/(8.)«B0J«A2P^5(ETA45<ETA4) 

CJ51P  =  B5/(B0)(A2P'ixA2P5(ETA'^X5<2XETA2) 

FUM1  =  3.5(THETA2-1. 

FUN2=1.-5.XTH£TA2 

SINEI2  =  SINEI5<X2 

A1  =  A2PXCJ25(FUN1 

A0=-A1/ETA3 

A2  =  3.3fA2PXCJ2XSINEI2 

FUN5  =  l.-ll.xTHETA2-(40.XTHETA'i)/FUN2 

FUN6=  -FUNl-(8.XTHETA'i)/FUN2 

FUNA=THETA2/SINEI2 

FUN22=FUN2XX2 

CJ21P2  =  CJ21PX3<2 

E01P  =  ((E2P5(ETA2)X(3.XCJ21P2XFUN5-10.X 
1CJ41PXFUN6))/(24.XCJ21P) 

E21P=-2.xE01P 

E31P  =  ((35.XCJ51P5(E2P25(ETA25<SINEI)XCFUN2-(16.XTHETA4)/FUN2))/(96.X 
1CJ21P) 

EllP  =  -.7  5D0xE31P+((  .2SD05(ETA25<SINEI)x(CJ31P+.3125D05tCJ51P 
1X(4.+3.»E2P2)X 
2(l.-9.xTHETA2-(2<4.XTHETA<+)/FUN2)))/CJ21P 

CI0=-(E2PXTHETA)/(ETA2XSINEI) 

CI2  =  CJ21P3<THETAJ(SINEI5«1.5D0 

CI 1  =  E2PXCI2X. 66666667  DO 

FUN7=(-.5DQXETA3XCJ21P)/E2P 

CL21P  =  (ETA3/CJ21P)?((  .  25D0XCJ21P2XFUN5- .83333333D0XCJ'ilPxFUN6  ) 

CL12P  =  CNOX(1.  +  1.5DOXCJ21PX£TAXFUN1+.0937  5D05(CJ21P2XETAX(-15.+I6. 
lxETA+25.xETA2+(30.-96.xETA-9  0.3(ETA2)XTHETA2+(105.+144.xETA+25. 
25(ETA2)3(THETAA)+.9  375D0XCJ41PXETAXE2P2X(3.-30.XTHETA2+35.XTHETA'+)) 

CL22P=.5D0XCN0XCN2 

G21P  =  (l./(2<+.5(CJ21P))x(-3.xCJ21P2x(2.+E2P2-ll  .  x(2  . +3  .  XE2P2)xTHETA2 
l-<40.X(2.+5.-XE2P2)XTHETA4/FUN2-40Q.5(E2P2XTHETA6/FUN22)  +  10.xCJ<4l?x 
2(2. 

3  +  E2P2-3.x(2.+5.;«E2P2)XTHETA2-8.X(2.+5.xE2P2)XTHETA4/FUN2-8  0.xE2P2^ 
4THETA6/FUN22)) 

G12P=Cn0X(-1.5D0xCJ21PXFUN2+.0937  5D0xCJ21P2X(-35.+24.xETA+25. 
1XETA2+ 

2(90.-192.?(ETA-126.xETA2)5(THETA2+(385.+36  0.xETA+'45.xETA2)5(THETA4)  + 
3.3125D05(CJ41P5((21.-9.XETA2+(-27  0.  +  126.5(ETA2)XTHETA2+(385.-189. 
4XETA2)?JTHETA4)) 

H2=1.5D0XCJ21PXTHETA 

u T  —  _2    3^H2 

H1=.66666667D0XE2PXH2 

H31P=((35.XCJ51PXE2P2XE2PXTHETA)/(144.XCJ21P))XC . 5D0/SINEIX ( FUN2-( 
116.XTHETA4)/FUN2)+SINEIX(5.+(32.XTHETA2)/FUN2+8  0.XTHETA4/FUN22)) 

H11P=-.25X 
1  H31P+i:C.25D0XE2P3<THETA)/(CJ21PXSINEI))x(CJ31P+.3125D0XCJ51Px 

2(4.+3.xE2P2)X(l.-9.XTHETA2-(24.XTHETA4)/FUN2)+1.875D0XCJ51P 
3XSINEI2X 
4(4.+3.xE2P2)x(3.+(16.xTHETA2)/FUN2+(40.xTHETA4)/FUN22)) 
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H21P  =  (E2P2XTHETA)/(12.3<CJ21P)X(-5.XCJ21P2X(11.  +  (3  0.XTHETA2)/FUN2+ 
l(200.XTHETA4)/FUN22)  +  10.XCJ'4lPX(3.  +  (16.XTHETA2VFUN2+(40.XTHETA4)/ 
2FUM22)) 

H12F=CN05(THETAX(-3.XCJ21P+.375D0XCJ21P2x(-5.  +  12.XETA+9.XETA2+ 
l(-55.- 
236.xETA-5.xETA2)xTHETA2)  +  1.25D0XCJ41P3<(S.-3.xETA2)X(3.-7.XTHETA2)) 

AID=CJ51P/CJ21P 

AID2=FUN2-(16.XTHETA^)/FUN2 

C1  =  35./38<4.5<AID5«ETA3XE2PXSINEIXAID2 

AID3=THETA2/SINEI 

AID4  =  THETA25(SINEI 

E2P3=E2P2XE2P 

C2  =  35./1152.XAIDX((-E2PXSINEI5((3.+2.XE2P2)  +  E2P3XAID3)XAID2+ 
12.5(E2P35<AID'+x(5.  +  (32.xTHETA2)/FUN2+(80.XTHETA4)/FUN22)) 

C3  =  l.-9.xTHETA2-(2'+.xTHETA4)/FUN2 

AID5=CJ31P/CJ21P 

C4=.25D0xAID5X(-E2PXAID3)  +  5./64.XAIDX(-E2P5(AID3x(4.+3.XE2P2)  + 
lE2PxSINElx(26.+9.xE2P2))xC3-15./32.XAIDxE2PXAID«tX(4.+3.xE2P2)X 
2C3.  +  (16.XTHETA2)/FUN2+(40.XTHETA'+)/FUN22) 

C5=E2P/(1 .+ETA3)X(5.-E2P2X(3.-E2P2)) 

C6  =  (E2Px(-32.+31.x(E2P2xE2P2)))/((4.+3.xE2P2)  +  ETAx(4.+9.3(E2P2)) 

C7=.25D0XAID5x3INElxC5+5./6a.xC3XAIDXETA2XSINElxC6 

C8  =  -.25D0xAID5xETA3XSINEI-5./(5'4.XAIDxETA3xSINElx(<4.+9.xE2P2)XC3 
0510  T=DT 

CL2P  =  CL12PxDT+CL22PxDTX5f2+CL02P  +  RNDXDTXDT 

CL2P=AM0D(CL2P, 6. 283 18530717  9600) 

IF(CL2P)520,530,530 
520  CL2P=CL2P+6.2831853071796D0 
530  G2P=G12PXDT+G02P 

H2P=H12PXDT+H02P 

SINEG=SIN(G2P) 

C0SING=COS(G2P) 

D1E=SINEGX(SINEGX(E31PXSINEG+E21P)+E11P)+E01P 

HlP=((H31PXSINEG+H21P;xSINEG+HllP)XC0SIMG+H2P 

GPLP=G2P+CL2P+.5D0X(CL21P+G21P)XSIN(2.XG2P)+(C1+C2)XC0S(3.X  G2P) 
1+(C4+C7)XC0SING 

CL1P=CL2P 

U=CL2P 
100  DELTAU=(U-E2PXSIN(U)-CL2P)/(1.-E2PXC0S(U)) 

U=U-DELTAU 

I F(ABS(DELTAU)-1. £-10)200, 100, 100 
200  U=U-(U-E2PXSIN(U)-CL2P)/(1.-E2PXC0S(U)) 

E  =  U 

SINE1P=SIN(E) 

C0SE1P=C0S(E) 

G1P=G2P 

ADIVR=1./(1.-E2PXC0SE1P) 

SINF1P=ADIVRXETAXSINE1P 

C0SrlP=ADIVRx(C0SElP-E2P) 

F1P=ATAN2(SINF1P,C0SF1P) 

I F(ABS(F1P-CL2P) -3. 141 5926 53589800)220,  21 0,210 
210    STOP 
220    F'jfl3  =  (l.+CN2XT)X3(.  66666667  DO 

C0SFG=C0S(2.X(G1P+F1P)) 

SINFG=SIN(2.X(G1P+F1P)) 

ADIVR2=ADIVR^X2 

ADIVR3  =  ADIVR>«>:3 

CI=CI2P+CIO^D1E+CI1XSINF1PXSINFG+(2.XCI1XCOSF1P+CI2)XCOSFG 

FUr;3  =  FlP-CLlP  +  E2PxSINFlP 
8  018    H  =  H1P+(2.^M1;-^'COSF1P+H2)XSINFG-H1XSINF1PXCOSFG+H3XFUN8 

KFUN=H/6.28313530717  96D0 

FUM9=KFUN 

H=H-FUN9X6 .28  318  530717  96  00 

IF(H)8022, 8023, 3023 

8022  H=H+6. 283185307179600 

8023  A=A2P/FUN3+A0+(A1+A2XC0SFG)XADIVR    +    AOTXOT 
AI06=A0IVR2XETA2+ADIVR 
AI07-SIN(2.XG2P+F1P) 
AID8=SIN(2.XG2P+3.XF1P) 
01=.25D0xCJ21Px(6.x(5.XTHETA2-l.)xFUN8+(3.-5.xTHETA2)X(3.xSINFG+ 

13.XE2PXAID7    +E2P?(AID8    )) 
D2=.2500XCJ21PX(2.X(3.XTHETA2-1.  )X(AID6  +  1  .  )xSINFlP+3  .  x(  1 . -THETA2)3< 
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1    <(-AID6+l. )XAID7+(AID6+.33333333D0)XAID8)) 

AID9=C0S(2.XG2P+F1P) 

AID10=COS(2.XG2P+3.XF1P) 

D3  =  -ETA2X.5DCXCJ21P5<(1.-THETA2)X(3.XAID9+AID10) 

ETA6I=1./(ETA3XETA3) 

D«t  =  ETA6I>5(C5  +C0SF1PX(3.+E2P5(C0SF1PX(3.+E2P5(C0SF1P))) 

D5  =  ETA6IX(E2P+C0SF1P5((3.+E2PXC0SF1P5((3.  +  E2PXC0SF1P))) 

D6=    ETA25(CJ2i<.5D0x((3.xTHETA2-l.))5D'i+3.x(l.-THETA2)5(D5XCOSFG)  +  D3 

GAL=GPLP+D1+(E2PXETA2)/(1.+ETA)XD2 

CE=(E2P-1.)X(1.+CN2XDT)5(3«.66666666666667D0  +  1.+D1E+D6    +    ESDxDT 

EDL=.5D05(E2P5<CL21PXSIN    (2.  XG2P)+CSXCOSING+E2P5(C1XCOS    (3.XG2P)- 
1      ETA5XD2 

AIDlif  =  SIN(CL2P) 

AID15=C0S(CL2P) 

ESL=CE  XAID14+EDLXAID15 

ECL=CE   XAID15-EDL5(AID14 

CE  =  SQRT(ECLJ<ECL  +  ESLXESL) 

CL=ATAN2(ESL,ECL) 

G=GAL-CL 

G=AMOD(G,6.28  318  53  071796D0) 

IF(G)302^,3025,S025 
S02<t  G  =  G  +  6.2831S53071796DO. 
8025  RETURN 

END 

X 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  XX  XX 

X  XX    SUBROUTINE  POSVEL     XX 

X  XX  XX  ' 

X  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX  I 

X  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

X 

X 

SUBROUTINE  POSVEL  (A,CE, CI , CL, G, H, BO, Rl, R2, R3, VI, V2, V3,R, VEL) 

CALL  NWTRPH  (CL,CE, E, SINEP,COSEP) 

HSIN=SIN(H) 

HCOS=COS(H)  I 

GSIN=SIN(G) 

GCOS=COS(G) 

CISIN=SIN(CI) 

CICOS=COS(CI) 

A1I=HC0SXGC0S-HSINXCIC0SXGSIN 

A12=-HC0SXGSIN-HSINXCIC0SXGC0S 

A21=HSINXGC0S+HC0SXCIC0SXGSIN 

A22=HC0SXCIC0SXGC0S-HSINXGSIN 

A31=CISINXGSIN 

A32=CISIMXGC0S 

FUN=SQRT(1.-CEXCE) 

Rl=AX(Allx(C03EP-CE)+A12x(FUNXSINEP)) 

R2=AX(A21X(C0SEP-CE)+A22X(FUNXSINEP)) 

R3  =  A5((A31XCC0SEP-CE)+A32X(FUNXSINEP)) 

R=Ax(l .-CE^CdSEP) 

F'JMl  =  (SQnT(BOxA))/R 

V1  =  FUN1);(A11X(-SINEP)  +  A12XFUN^C0SEP) 

V2  =  FUMl5«(A2lx(-SINEP)+A22?(FUM5(COScP) 

V3  =  FUtllx(A3l5<(-SinE?)+A32XFUNXC0SEP) 

VEL=FUNlx3qRT(l.-(CExC0S£P)xx2J 

RETURN 

END 
X 
X 

X  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

X  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

X  XX  XX 

X  XX    SUBROUTINE  NWTRPH     xx. 

X  XX  XX 

X  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

X  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 

X 
X 

144 


SUBROUTINE  NWTRPH  (QL , B, E2, SE, CE) 

NEWTON  RAPHSON  SOLUTION  TO  KEPLERS  EQUATION  SUBROUTINE 

E1  =  QL  +  (B5<SINCQL))/(1.-BXC0S(QL)) 

SE=SIN(E1) 

CE=C0S(E1) 

E2=E1+(QL+BXSE-E1)/(1.-BXCE) 

IF(ABS(E2-E1)-  l.E-8  )'i,4,3 

E1  =  E2 

GO  TO  2 

SE=SINCE2) 

CE=C0S(E2) 

RETURN 

END 


XX  XX 

XX        FUNCTION  XSPA       XX 

XX  XX 

XXXXXHXXXXXXXXXXXXXXXXXXXXXXXX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


THE  FOLLOWING  FUNCTIONS  DETERMINE  THE  FIRST  ORDER  SHORT 
PERIODIC  VARIATIONS  AS  A  FUNCTION  OF  THE  MEAN  ELEMENTS. 


FUNCTIO 
DIMENSI 
REAL  J( 
DATA  J, 
F  =  XM( 

1+(1.25D 

2XM( 

26))+(13 
P  =  XM( 
RR  =  P/ 
XSPA=  J 

HXM(3)) 

2(1.-3./ 
RETURN 
END 


N  XSPA(XM) 

ON  XM(6) 

4) 

R/1.,0.00108276DO,-0.00000255DO,-0.00000156DO,6378.165DO/ 

6)+(2.xXM(2)-XM(2)XX3./4+5./96.xXM(2)xx5. )XSIN(XM(6)) 

0XXM(2)XX2.-11./24.XXM(2)XX4.+17./192.XXM(2)XX6.;XSIN(2.X 

./12.xXM(2)xx2.-43./192.xXM(2)XX5.)XSIN(3.XXM(6)) 

1)X(1.-XM(2)XX2.) 

(1.+XM(2)XC0S(F)) 

(2)XRXR/XM(1)X((XM(1)/    RR      )XX3  .  X( (1 . -3 ./2 . XSIN(XM( 3) )XSIN 

)  +  3./2.xSIN(XM(3))xSIN(XM(3))xC0S(2.x(XM(<4)  +  F)))- 

2.XSIN(XM(3))XSIN(XM(3)))X(1.-XM(2)XX2.)XX(-1.5D0)) 


xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

XX  XX 

XX        FUNCTION  XSPE       XX 

XX  XX 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 


FUNCTION  XSPE(XM) 

DIMENSION  Xri(6) 

REAL  J (4) 

DATA  J,R/1  .,0.001 08276 DO, -0.00 00025500,-0. 00000156 00,6  378.16500/ 

F  =  XM(6)  +  (2.xXM(2)-XI1(2)xx3./4  +  5./96.xXN(2)xx5.  )XSIN(XM(6)  ) 
l  +  (1.25D0xXM(2)xx2.-ll./2<+.xxM(2)xx4.  +  17  ./192  .  XXM(2)XX6  .  )XSIN(2.x 
2XM( 
36))+(13./12.xXM(2)xx2.-43./192.xXM(2)xx5.)xSIN(3.xXM(6)) 

P  =  XM(l)x(l.-XM(2)xx2.  ) 

XSPE  =  J(2)/2.X(R/P)XX2.X(1.-1.5D0XSIN(XM(3))XSIN(XM(3))) 
1X((1.+1.5D0XXM(2)XXM(2)-(1.-XM(2)XXM(2))XX1.5D0)/XM(2)+3.X(1.+ 
2XM(2) 

3xxrU2)/4.)xCOS(F)  +  1.5D0XXM(2)xCOS(2.xF)  +  XM(2)XXM(2)/<4.XCOS(3.XF)) 
<++3./S.xj(2)x(R/P)xx2.xSIN(XM(3))xSIN(XM(3))X((l.+ll./4.xXM(2)xx2. ) 

5  XC0S(2.XXM(<4)+F)+XM(2)XXM(2)/4.XC0S(2.XXM(4)-F)  +  5.XXM(2)X 

6  C03(2.xXM(4)+FX2.)+(7.+17./4.xXM(2)xXMC2))/3.xC0S(2.XXM(4)+3.XF) 
7+1.5D0XXM(23XC0S(2.xXM(4)+4.XF)+XM(2)xXM(2)/4.xCOS(2.xXM(4)+5.xF) 
8+1.5D0xXM(2)xcOS(2.xXM(4))) 
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RETURN 

END 
X 
X 

X  XX  XX 

X  XX        FUNCTION  XSPI       XX 

X  XX  XX 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 
X 

FUNCTION  XSPI(XM) 

DIMENSION  XM(6) 

REAL  J(«t) 

DATA  J, R/1., 0.00 108276 DO, -0.00 0  0025500,-0. 0  0 000 156 DO, 6 378. 16 5D0/ 

F  =  XM(6)  +  (2.xXM(2)-XM(2)xx3./<*+5./96.xXM(2)Xx5.  )XSIN(XM(6)) 
1+(1.25D0XXM(2)XX2.-11./24.XXM(2)XX4.+17./192.XXM(2)XX6 .)XSIN(2.X 
2XM( 
36))  +  (13./12.xXM(2)xx2.-<43./192.xXM(2)xx5.)xSIN(3.xXM(6)) 

P  =  XM(1)X(1.-XM(2)XX2.) 

XSPI  =  3./8.xj(2JX(R/P)xx2.XSIN(2.5«XM(3))X(XM(2)xC0S(2.XXM(4)  +  F) 
1  +C0S(2.X(XM(4)  +  F))+XM(2V3.XCQS(2.XXM(4)  +  3.XF)) 

RETURN 

END 
X 
X 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  XX  XX 

X  XX       FUNCTION  XSPW      xx 

X  XX  XX 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 
X 

FUNCTION  XSPW(XM) 

DIMENSION  XM(6) 

REAL  J(^) 

DATA  J, R/ 1., 0.0 01 08276 DO, -0. 00 000255D0, -0.0 0  00 0156 00,6  378.16500/ 

F  =  XM(6)+(2.xXM(2J-XM(2)XX3./^+5./96.XXM(2)xx5.)XSIN(XM(6)) 
1+(1.25D0XXM(2)XX2.-11./2^.XXM(2)XX4.+17./192.XXM(2)XX6.)XSIN(2.X 
2XM( 
26))+(13./12.xxM(2)xx2.-^3./192.xXM(2)xx5.)xSIN(3.xXM(6)) 

P  =  XM(1)X(1.-XM(2)X5(2.) 

XSPW  =  3./4.xj(2)X(R/P)xx2.x(«i.-5.xSIN(XM(3))XSIN<XM(3))) 
1X(F-XM(6)+XM(2)XSIN(F))+1.5D0XJ(2)X(R/P)XX2.X(1.-1.5D0XSIN(XM(3)) 
2XX2.) 

3X((l.-XM(2)XXM(2)/'^.)/XM(2)XSINCF)+.5D0xSIN(2.xF)+XM(2)xSIN(3.XF) 
4/12 
5.)-1.5D0xj(2)x(R/P)XX2.x((SIN(XM(3))xSIN<XM(3))/<i.+XM(2)xx2./2.x( 

61.-15./8.XSIN(XM(3))XX2.  )  )/XM(2)xSIN(  2  .  XXM( ';)  +  F)  +  XM(  2)/16  .  XSIN 
7(XM(3))XX2.XSIN(2.XXM(4)-F)+.5D0X(1.-2.5D0XSIN(XM(3))XX2..)XSIN(2.X 

8(XM(4)  +  F))-C7./12.xSIN(XM(3))xx2.-XM(2)XXM(2)/6.x(l.-19./8.)f5IN 
9(XM(3))XH2.  ))/XM(2)xSIN(2.xXM(4)  +  3.?<F)-3./3.  XSIN(XM(3))XX2.XSIN 
11(2.;5XM(4)+4.XF)-X:i(2)/16.XSIN(XM(3))XX2.XSIN(2.XXM(4)+5.XF)) 
12-9./16.xj(2)X(R/P)XX2.xsiN(XM(3))xx2.xSIN(2.xXM('4)) 
RETURN 
END 
X 
X 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  XX  XX 

X  XX       FUNCTION  XSPO      XX 

X  XX  XX 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 
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FUNCTION  XSPO(XM) 

DIMENSION  XM(6) 

REAL  J (4) 

DATA  J, R/1., 0.001 08276 DO, -0.00000255D0, -0.00000156 DO, 6 37S.I65D0/ 

F  =  XM(6)  +  (2.5<XM(2)-XM(2)X3<3./A+5./96.5(XM(2)5<5(5.  )XSIN(XM(6)) 
l  +  (1.25D0XXM(2)XX2.-H./2'i.XXM(2)x»4.+17./192.XXM(2)xx6.)xSIN(2.x 
2XM( 
36))  +  (13./12.xXM(2)5(5«2.-43./192.xXM(2)XX5.)xSIN(5.xXM(6)) 

P  =  XM(1)^(1.-XM(2)X)(2.) 

XSPO  =  -1.5DOXJ(2)3((R/P)X5<2.XCOS(XM(3))5((F-XM(6)+XM(2)XSIN(F)-XM(2) 
l/2.5tSIN(2.xXM(a;  +  F)-SIN(2.x(XM(4)  +  F))/2.-XM(2)/6.xSIN(2.xXM(4)  +  3.x 
2F)) 

RETURN 

END 
X 
X 
X  3(9(X)(X3(XXKXXXXXXXXXXXXXX9(HXXXXX 

X  )(X  XX 

X  XX        FUNCTION  XSPM       XX 

X  XX  XX 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 
X 

FUNCTION  XSPM(XM) 

DIMENSION  XM(6)- 

REAL  JCi) 

DATA  J, R/1., 0.00 108276 DO, -0.0000025500,-0. 000 0  0156 DO, 6 378. 165D0/ 

F  =  XM(6)+(2.xXM(2)-XM(2)XX3./4+5./96.xXM(2)xx5.)xSIN(XM(6)) 
l  +  (1.25D0XXM(2)xx2.-ll./2<+.XXM(2)xx<4.  +  17./192.xXM(2)XX6.)xSIN(2.X 
2XM( 
36))  +  (13./12.xXM(2)xx2.-<i3./192.xXM(2)xx5.)xSIN(3.xXM(6)) 

P  =  XM(1)X(1.-XM(2)XX2.) 

XSPM  =  -1.5D0XJ(2JX(R/P)XX2.XSQRT(1.-XM(2)XX2.)/XM(2)X((1.-1.5D0X 
1SIN( 

2XM(3))XX2.)X((1.-.25D0XXM(2)XXM(2))XSIN(F)+XM(2)/2.XSIN(2.XF)+ 
3XM(2) 

4XXM(2)/12.XSIN(3.XF))+.5D0XSIN(XM(3))XSIN(XM(3))X((-.5D0)X(1.+ 
51.25D0 

6xXM(2)xXM(2))xSIN(2.xXM(4)  +  F)-XM(2)xXM(2)/8.xSIN(2.xXM('i)-F)  +  7./6. 
7x(l.-XM(2)xx2./28.)xSIN(2.xXM(<i)  +  3.xF)+.7  5D0XXM(2)xSIN(2.xXM(4)+<i. 
8  xF 

9)  +  XM(2)xx2./8.xSIN(2.xXM('4)  +  5.XF)))  +  9./16.xj(2)X(R/P)xx2.x 
1SQRT(1.-XM(2)XX2.)XSIN(XM(3))XSIN(XM(3))XSIN(2.XXM(4)) 

RETURN 

END 
X 
X 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  XX  XX 

X  XX        FUNCTION  XSPR       XX 

X  XX  XX 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

X 

FUNCTION  XSPR(XM) 

DIMENSION  XM(6) 

REAL  J(A; 

DATA  J,R/1.,0.00108276DO,-0.00000255DO,-0.0  0  00  0156DO,6  37  8.165DO/ 

F  =  XM(6)  +  (2.xXM(2)-XM(2)XX3./«4+5./96.xXM(2)xx5.  )XSIN(XM(6)) 
1+(1.25D0XXM(2)XX2.-11./24.XXM(2)XXA.+17./192.XXM(2)XX6. )XSIN(2.X 
2XM(6n  +  (13./12.xXM(2)XX2.-<43./192.xXM(2)XX5.)x3IN(3.xXM(6)) 

P  =  XM(1)X(1.-XM(2)X3<2.  ) 

RR  =  P/(1.+XM(2)XC0S(F)) 

XSPR  =  -.5D0XJ(2)X(RXX2./P)X(1.-1.5D0XSIN(XM(3))XSIN(XM(3)))X(1.+ 
1XM(2) 

2/(l.+SQRT(l.-XM(2)xx2.))xC0S(F)+2./(SQRT(l.-XM(2)xx2.))xRR/XM(l 
3))+.25XJ(2)xRxx2./PxSIN(XM(3))xx2.xC0S(2.xXM(<4)  +  2.XF) 
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RETURN 

END 
X 
X 

X  XX  XX 

X  XX      SUBROUTINE  MEAN       XX 

X  XX  XX 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 
X 


SUBROUTINE  MEANCXOSCXM) 

COMPUTES  THE  MEAN  ORBITAL  ELEMENTS  FROM  THE  OSCILATORY  ELEMENTS 
AND  THE  FIRST  ORDER  PERIODIC  VARIATIONS  BY  THE  RELATIONS, 

XM  =  XOSC  -XSP(XM) . 

DIMENSION  X0SC(6),XM(6),XMN(6),XSP(6),T0L(6),DELXM(6) 

DATA  TOL/.OIDO, .00001  DO,. OOOOI DO,. 00001  DO,. 00001  DO,. 00001  DO/ 

KOUNT  =  0 
X...   INITIAL  GUESS  FOR  MEAN  ELEMENTA. 

DO  10  1=1,6 
10  XM(I)  =  XOSC(I) 
X...   CALCULATE  THE  SHORT  PERIODIC  VARIATIONS. 
20  XSPCl)  =  XSPA(XM)- 

XSP(2)  =  XSPECXM) 

XSP(3)  =  XSPKXM) 

XSP(4)  =  XSPW(XM) 

XSP(5)  =  XSPO(XM) 

XSP(6)  =  XSPM(XM) 

KOUNT  =  KOUNT+1 

DO  50  J=l,6 

XMN(J)  =  XOSC(J)-XSP(J) 

DELXM(J)  =  ABS(XMN(J)-XM(J)) 

IF(DELXM(J).GT.TOL(J))  XMCJ)  =  XMN(J) 
50  CONTINUE 

DO  60  K=l,6 

IF(DELXM(K).GT.TOL(K).AND.-KOUNT.LT.30)  GO  TO  20 
60  CONTINUE 

70  IF(KOUNT.GT.30)  WRITE(6,100)  DELXM 
90  CONTINUE 
100  F0RMAT(/T5, 'CONVERGENCE  OF  ONE  OR  MORE  MEAN  ELEMENTS  WAS  NOT  MET.' 
l/'THE  ABSOLUTE  DIFFERENCE  IN  THE  ELEMENTS  ARE  SHOWN. '//7E12 . 5) 

RETURN 

END 
X 
X 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  XX  XX 

X  XX     SUBROUTINE  GEOP      XX 

X  XX  nr. 

X  XXXXXXXXXXXXXXXXXXXXXXXXXXXXXK 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 
X 

SUBROUTINE   GEOP    (XM,DXDT) 
DIMENSION    XM(6)    ,    DXDT(6) 
REAL    N,NCONS,NCONST,NCON 
REAL    JC^) 

DATA    J, R/1., 0.00 108276 DO, -0. 00 0  00255DO, -0.00000156 DO, 6 378. 165D0/ 
DATA    U/398600./ 
N    =    SQRT(U/(XM(1)XX3.)) 
P    =    XM(1)X(1.-XM(2)XX2.) 
DXDT(l)    =    0. 

S3    =    SIN(XM(3))XSIN(XM(3)) 
XS    =    l.-XM(2)XX2. 
CS    =    C0S(XM(3)) 

DED    =         (-3.  )/32.XNXj(2)XJ(2)x(R/P)xx4.xS3X(l<+.-15.xS3)xXM(2)X 
1XSXSIN(2.XXMC^)) 
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DEDO    =    3./8.XNXJ(3)X(R/P)XX3.xSIN(XM(3))X(4.-5.xS3)XXSXC0S(XM(<4)) 

DEDOT    =    15./32.XNXJ((+)3((R/P)xx<+.XS3X(6.-7.3<S3)X(XM(2)5«XS) 
1XSIN(2.5(XM(<+)) 

DXDT(2)    =    DED-DEDO-DEDOT 

DID   =    3./6'4.5(NXJ(2)5(J(2)x(R/P)X5(4.3<SIN(XM(3)x2.)x(l'i.-15.xS3)x 
IXr-K  2 )  XX2  .  HSI N ( 2  .  stXi-K  «♦ ) ) 

DIDO    =    3./S.XM3<J(3)X(R/P)XX3.xC0S(XM(3))3<(4.-5.xS3)5«XM(2)3«C0S(XM(4 
D) 

DIDOT    =    15./64.XNXJ(4)X(R/P)X5«4.XSIN(2.XXM(3))X(6.-7.XS3)XXM(2)X 
1XM(2)XSIN(2.5«XM('^)) 

DXDTC3)    =    DID+DIDO+DIDOT 

DND    =    .75D0XNXJ(2)5«(R/P)XX2.3<('^.-5.XS3) 

DH    =    2.X(14.-15.5<S3)XS3-(2S.-158.XS3+135.5(S3XX2.)XXM(2)X3«2. 

DWDO  =  3./16.XN>(J(2)5(^2.5<(R/P)xx<i.x(48-103.xS3  +  215  ./"i  .  XS3X5(2  .  + 
l(7.-A.5D0XS3-'i5./S.xS3X5(2.)*XM(2)xx2.+6.x{l.-1.5D0XS3)x(4.-5.xS3)x 
2     SQRT(XS)-.25DOXDW3(COS(2.»XM(|^))) 
»...   PRINT  *,    XM(2),  XMC3),  Xf-K-^) ,  R,  P,  N,  J(  3) ,  S3 

DWDOT  =  3./8.)(NXJ(3)X(R/P)3<X3x((4.-5.xS3)x(S3-XM(2)5(3f25(COS(XM(3) 
l)X5(2)/(Xf-1(2;xSIH(XM(3)))  +  2.xsiN(XM(3)))«(13.-15.xS3)XXM(2nxSIN(XM 

DHW  =  S3>((6.-7.xS3)-(6.-35.xS3  +  31.5xS3xx2.)XXM(2)5(X2. 

DllDOTT  =  15./32.XNXJ(<i)J((R/-P)XK4.x(16.-6  2.5<S3+<i9.xS3xx2.  +  .75D05((2'+.- 
18'4.XS3+63.XS3XX2.)5«XM(2)X5(2.+DHW5<C0S(2.XXM(4))) 

DXDT(4)  =  DWD+DWDO+DWDCT-DWDCTT 

DOD  =(-1  .5D0)XNXJ(2)X(R/P)XX2.XCS 

DODO  =  1.5D0)(NXJ(2)X;(2.)((R/P)X)(<*.XCSX(2.25D0  +  1.5XSQRT(XS)-S3 
l3f(2.5D0+2.25D0xSQRT(  XS)  )+XM(2)x>C2  ./4  .  *(  1  . +1 .  25D0XS3)-XM(2)3(X2  . 
2/8.x(7.-15.xS3)xCCS(2.5(XM('4n) 

DODOT  =  3./8.XN5(J(3)x(R/P)xx3.X(15.XS3-4)XXM(2)XCS/SIN(XM(3))5(SIN( 
1XM(<4)) 

DODOTT  =  15./16.XNXJ(<i)x(R/P)X3f<i.xCSx((<i.-7.»S3)X(l.+1.5D0XXM(2)xx2. 
2)-(3.-7.xS3)xxM(2)xx2.xC0S(XM('4)x2.)) 

DXDT(5)  =  DOD-DODO-DODOT+DODOTT 

NCON  =  3./8.xNXj(2)XX2. 

NCOMS  =  NXJ(<;)!<(R/P)xx<4. 

NCONST  =  3./8.5(Nxj(3)X(R/P)XX3. 

DM  =  NX(1.+1.5D0XJ(2)X(R/P)XX2.X(1.-1.5D0XS3)XSQRT(XS)) 

DMD=I.5D0XNXJ(2)XX2?<(R/P)xx<+X((l.-1.5D0XS3)XX2XXS+(1.25D0X(1.- 
12.5D0X 

lS3  +  13./8.xS3XS3)  +  5./8.x(l.-S3+5./8.xS3xS3)xXM(2)xx2+S3/16.x(l<i.- 
2     15.xS3)x(l.-2.5D0xxri(2)xx2)xCOS(2.xXM(i;)))xSQRT(XS)) 

DM1 =(3. -7 .5D0XS3+47./8.xS3xS3+(1.5D0-5.xS3+117./6.xS3xS3)xXM(2)XX2 
l-(l.+5.xS3-101./8.xS3xS3)/8.xXM(2)xx<i.) 

DM2  =  XM(2)X3(2./8.xS3x(7  0.-123.xS3+(56  . -66  .  XS3)XXM(2)XX2.  )XC0S(2  .  X 
lXM(4))  +  27./128.xXM(2)xx<i.xS3xS3xC0S(4.xXM(<4)) 

DMDO  =  NC0NX(R/P)XX4./(SQRT(XS))X(3.XDM1+DM2) 

DMDOT  =  NC0NSTXSIN(XM(3))X(4.-5.XS3)X(1.-4.XXM(2)XX2.)/XM(2)X 
1SQRT(XS)X 
lSIN(XM(4))-<+5./128.xNCONSx(8.-40.xS3+35.xS3xS3)xXM(2)xx2.xSQRT(XS) 

DMD0TT  =  15./6«t.xNC0HSxS5x(6.-7.xS3)x(2.-5.xXM(2)X5<2.  )XSQRT(XS)XC0S( 

12.XXM(4)) 
DXD~(6)    =    DM+DMD+DMDC-DMDOT+DMDOTT 

RETURN 

e:id 

X  • 

X  XXXXXXXXXXXXXXXX5€3tXX3<XXXXXXXXX 

»  )-:X3<XXXXXXXXXXXXXXXXXX3(XXXXXXXX 

X  XX  XX 

X  XX      SUBROUTINE  SOLOC      xx' 

X  XX  XX 

X  XXXXXXXXX3(5(XXXXXXXXXXXXXXX5<X=<X 

X  XXXXXXXXXXXX3(X^XXXXXXXXXXXXXXX 

X 
X 

SUBROUTINE   SOLOC    (    TVE    ,    T    ,    BLA    ,    DS    ,    AS    ,    AB    ) 
DATA    EPSl/    23.5D0   /    ,    PI   /    3.1<+159D0/,    TOL/    O.OOOIDO   / 
PIl    =    0.5D0XPI 
PI2    =    2.    X    PI 
PI3    =    1.5D0X   PI 
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RPD  =  PI  /  180. ODO 
DELT  =  T  -  TVE 
EPS  =  EPSl  X  RPD 
CLON  =  0.985647D0  x  RPD  X  DELT 
IF  (CLON.LT.  0.0  )  CLON  =  PI2  +  CLON 
IF(CL0N.GT.PI2)  CLON  =  AM0D(CL0N,PI2) 
SDS  =  SINCEPS)  X  SIN  (  CLON  ) 
DS  =  ASIN  (  SDS  ) 
BLAR  =  BLA  x  RPD 
DO  10    1=1,5 
J  =  I  -  1 

ARC  =  0.5D0  X  J  X  PI 
TST  =  ARG  -  CLON 
TSTl  =  ABS(  TST  ) 
IF  (  TSTl  .  LE  .  TOL  )  GO  TO  20 
10  CONTINUE 

ETAl  =  TAN(DS)  /  TAN(EPS) 

ETA  =  ASIN(ETAl) 

ETA  =  ABS(ETA) 

IF  ((  CLON  .  LT  .  PIl  )  .  AND  .  (  CLON  .  GT  .  0 . 0) )  AS  =  ETA 

IF(CCLON.LT.PI) .AND.(CLON.GT.PIl))  AS  =  PI  -  ETA 

IF  ((CLON. LT. PIS) .AND. (CLON. GT. PI))  AS  =  PI  +  ETA 

IF  ((CLON. LT.PI2). AND. (CLON. GT.PI3))  AS  =  PI2  -  ETA 

GO  TO  21 

20  AS  =  CLON 

21  A3  =  AS  +  BLAR 
RETURN 

END 
X 
X 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  XX  XX 

X  XX    SUBROUTINE  DATPRP     XX 

X  XX  XX 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 
X 

SUBROUTINE  DATPRP  (  X,Y,XX,YY  ) 

COMMON  /  PREP  /  ASUN,  DSUN,  RAS,  DECS 

X  =  ASUN 

Y  =  DSUN 

XX  =  RAS 

YY  =  DECS 
.  RETURN 

END 
X 
X 

X  XXXXXXXXXXJJXXXXXXXXXXXXXXXXXXX 

X  XXXXXXXXXXXXJtXXXXXXXXXXXXXXXXX 

X  XX  XX 

X  XX     SUBROUTINE  PERIOD      XX 

X  XX  XX 

X  xxxxxxxxxxxxxxxxxxxxx;«xxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 
X 

SUBROUTINE  PERIOD  (  AO,  ES,  XI,  PA  ) 

REAL  J2 

COMMON  /  XXXX  /  lET 

DATA  J2  /  0.00108276DO/,  R/  6378.165D0/  ,  U/398600 ./, PI/3 . iai59D0/ 

PA  =  2.XPIX(SQRT(A0XA0XA0/U))X(1.»-(3.XJ2XRXRX(3.XSIN(XI)XSIN(XI) 
1    -2.0)/(4.XA0xA0x((l.-ESXES)xxl.5D0)))) 

IF  (  lET  .  EQ.  0  )   RETURN 

WRITE  (6,100)  PA 
100  FCRMAT(/10X, 'INITIAL  ANOMALISTIC  PERIOD=  ',F12.3,'  SEC',//) 

RETURN 

END 
X 
X 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
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X  XX  XX 

X  XX  SUBROUTINE   SALT  XX 

X  XX  XX 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 
X 

SUBROUTINE  SALT  (D,XM, RM,Z, RE  ) 

DIMENSION  XM(6) 

DATA  R0£/6  373.16  5DO/,EOE2/6.6  9342E-03/ 

RE  =  R0E/SgRT(l.+(E0E2/(l.-E0E2))xSIN(D)xSIN(D)) 

Z  =  RM  +  XSPR(XM)  -  RE 

RETURN 

END 
X 
X 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  XX  XX 

X  XX     SUBROUTINE  SATLOC      xx 

X  XX  XX 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 
X 

SUBROUTINE  SATLOC  (MTYPE, XM, DELT, DSUN, ASUN, DECS, RAS,GMGLT, R) 

DIMENSION  XM(6) 

DATA    B0/3986  0  3.254DO/,W/.7292115E-0^/,PI/3.1<+159DO/ 

AM    =    0.0 

PI02  =  0.5D0XPI 

PI2  =  2.XPI 

P2  =  PI2 

PI03  =  3.XPI/2. 

AG  =  AMOD  (  WXDELT  ,  PI2  ) 

A  =  XM(1) 

CE  =  XM(2) 

CI  =  XM(3) 

G  =  AMOD  (  XMCi)  ,  PI2  ) 

H  =  AMOD  (  XM(5)  ,  PI2  ) 

CALL  POSVEL(A,CE,CI,AM,G,H,B0,Rl,R2,R3,Vl,V2,V3,R,V) 

DECS  =  ATAN  (  R3/SQRT(  RlxRl  +  R2XR2  )) 

RAS  =  ATAN2  (  R2  ,  Rl  ) 

IF  •(  MTYPE  .EQ.  0  )  GO  TO  10 

ALP  =  TAN(DSUN)  X  COSCCI)  /  SIN(CI) 

ALP  =  ASIN(ALP) 

ALP  =  ABS(ALP) 

DECS  =  DSUN 

IF  (  MTYPE  .  EQ  .  1  )  DECS  =  -  DSUN 
X      Ul  =  ASIN(  SIN(DSUN)  /  SIN(XM(3))) 
X      U2  =  ASIN  (  SIN(DECS)/SIN(XM(3))) 

Tl  =  ASUN  +  PI02 

T2  =  ASUN  -  PI02 

WN  =  G 

DMA  =  H 

IF  (T2  .LT.  0.0  )  GO  TO  20 
■   IF  (  Tl  .GT.  PI2  )  GO  TO  30 

IF  ((OMA  .LE.  Tl)  .AND.  (CMA  .GE.  T2))  GO  TO  ^0 

GO  TO  50 
20  T3  =  PI2  +  T2 

IF  (((OriA.LE.Tl)  .AND.  (  OMA.  GE.  0.0))  .  OR  .  ( ( OMA  .  LE .  P2)  .  AND.  (OMA  .  GE.  T3 
1    )))  GO  TO  40 

GO  TO  50 
30  T3  =  Tl  -  PI2 

IF(((OMA.GE.T2).AND.(OMA.LE.PI2)).OR.((OMA.GE.0.).AND.(OMA.LE.T3 
1    )))  GO  TO  40 

GO  TO  50 

40  IF  (  CI  .  GE  .  PI02  )  GO  TO  41 
RAS  =  OMA  +  ALP 

IF  (  DECS  .  LT  .  0.0  )  RAS  =  OMA  -  ALP 
GO  TO  42 

41  RAS  =  OMA  -  ALP 
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IF  (  DECS  .  LT  .  0.0  )  RAS  =  OMA  +  ALP 

42 

IF  (  MTYPE  .  EQ  .  1  )  RAS  =  RAS  +  PI 
RAS  =  AMOD  (  RAS  ,  PI2  ) 
GO  TO  10 

50 

IF  (  CI  .  GE  .  PI02  )  GO  TO  51 

RAS  =  OMA  +  PI  -  ALP 

IF  (DECS  .  LT  .  0.0  )  RAS  =  OMA  +  PI  +  ALP 

GO  TO  42 

51 

RAS  =  OMA  +  PI  +  ALP 

IF  (  DECS  .  LT  .  0.0  )  RAS  =  OMA  +  PI  -  ALP 

GO  TO  42 

10 

ELNG  =  RAS  -  AG 

IF  (ELNG  .LT.  0.0  )  ELNG  =  PI2  +  ELNG 
ELNG  =  AM0D(ELNG,PI2) 

GMGLT  =  ASIN  (  0 . 9792D0XSIN(DECS)+0 .2028D0XCOS(DECS)XCOS( ELNG  - 
I    (291.XPI/180.))) 
RETURN 

X 
X 

END 

X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 

XX                                  XX 

X 

XX      SUBROUTINE  DRAG       xx 

X 

XX                                  XX 

X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 
X 
X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

SUBROUTINE  DRAG(  XM,  DEL,  TT,  DELMG  ) 

REAL  N, 10, II, 12, 13, 14, 15, 16, 17 

COMMON  /  SOLAR  /  F107 , FBAR, AKP, TVE 

COMMON  /  PREP  /  ASUN,DSUN, RAS, DECS 

DIMENSION  XM(6),  SUN(4)  ,  SAT(2)  ,  GE0(3)  ,  BF(8) 

DATA    BO/3986  03.254DO/,BLA/30./,PI/3.14159D0/,EOE/.00  335D0/ 

DATA    IYR/1981D0/,TVD/1./,THR/17./,TM/14./,TS/56.243D0/ 

PI2  =  2.  X  PI 

SUN(3)  =  23.5D0X(PI/180.) 

SUN(4)  =  1.0 

GEO(l)  =  F107 

GE0(2)  =  FBAR 

GE0(3)  =  AKP 

H  =  SQRT(B0/(XM(1)XX3)) 

TCOR  =  AMOD  (  XM(6)  ,  PI2  ) 

TCOR  =  TCOR  /  N 

TTT  =  TT  -  (  TCOR  /  86400.  ) 

TTR  =  1721044.  +  367XIYR  -  (7XIYR/4)  +  TVD  +  ((  THRX3600.  +  TMX60 
1     +  TS  )  /  86400.  )  -  2400001. ODO 

DELT  =  TTT  -  TTR 

CALL  SOLOC(TVE,TTT,BLA,DSUN,ASUN,ABLG) 

MTYPE  =  0 

CALL  SATLOC( MTYPE, XM, DELT, DSUN,ASUN, DECS, RAS, GMGLT, RM) 

STR  =  XM(6) 

X'<1(6)  =  0.0  0 

CALL  SALT(  DECS, XM, RM, HP,RE  ) 

XM(6)  =  STR 
X      CALL  DATPRP(  SUN( 1) , SUN(2) , SAT( 1) , SAT( 2)  ) 
X      CALL  ATMDEN(TTT, SUN, SAT, GMGLT, HP, GEO, SHO,DRDH,SH) 
X      RHOP  =  RHO  X  (  10.  XX9.  ) 
X      SHP  =  SH  X  O.OOIDO 

CALL  JAC6  0(DECS,DSUN,RAS,ASUN,BLA,HP,F107,RHOP,SHP) 

MTYPE  =  1 

CALL  SATLOC(MTYPE,XM, DELT, DSUN,ASUN, DECS, RAS, GMGLT, RM) 
X      CALL  DATPRP  (  SUN( 1) , SUN(2) , SAT( 1) , SAT( 2)  ) 
X  •    CALL  ATMDEN(TTT, SUN, SAT, GMGLT, HP, GEO, RHO, DRDH,SH) 
X      RHOMIN  =  RHO  X  (  10.XX9  ) 
X      SHMAX  =  SH  X  O.OOIDO 

CALL  JAC6 0( DECS, DSUN, RAS, ASUN,BLA, HP, Fl 07, RHOMIN, SHMAX) 

MTYPE  =  2 

CALL  SATLOC( MTYPE, XM, DELT, DSUN, ASUN, DECS, RAS, GMGLT, RM) 
X      CALL  DATPRP  (  SUN(l),  5UN(2),  SAT(l),  SAT(2)  ) 
X      CALL  ATMDEN(TTT, SUN, SAT, GMGLT, HP, GEO, RHO, DRDH,SH  ) 
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K      RHOMAX=RHO  X  (10.3<X9  ) 
X      SHMIN  =  SHXO.OOIDO 

CALL  JAC6  0(DECS,DSUN,RAS,ASUN,BLA,HP,F107,RHOMAX,SHMIN) 

A  =  SIN(DSUN)3«SIN(XM(3))XSIN(XM(<4))+C0S(DECS)X(C0S(XM(5)-ABLG)X 
1    C0S(XM(4))-CaS(XM(3))J(SIH(XM(5)-ABLG)^SIN(XM(':tn) 

B  =  SIN(DSUN)XSIN(XM(3))XC0S(XMC«4))-C0S(DECS)X(C0SCXM(5)-ABLG)3( 
1    SIN(XM(4))+C0S(Xr'1(3))»SIN(XM(5)-ABLG)XC0S(XM(^))) 

F  =  (RHOMAX  -  RHOMIN)  /  (RHOMAX  +  RHOMIN  ) 

RHOO  =  RHOP  -  0.5DOX(RHOMAX  +  RHOMIN)  X  F  X  A 

Z  =  XM(1)  X  XM(2)  /  SHP 

Zl  =  Zx(Z  -!.)-(  FXA  /  (1.  +  F  X  A  )) 

Z2  =  (  1.  /  (  1.  +  F  X  A  ))  -  Z  X  (  1.  -  0.5D0XZ)  -  0.5D0 

TMl  =  -ZlXZl  -  2.  X  ZXZXZ2 

CTHAVG  =0.0 

IF  (  TMl  .GT.  0.00  )  CTHAVG=(Z1+SQRT(TM1))/(ZXZ) 

SHNM  =  ABS(  SHMAX  -  SHMIN  ) 

HAVG  =  0.5D0X(SHMAX+SHMIN)+  (  SHNM  /  (SHMAX+SHMIN) )XAXCTHAVG 

BETA  =  1.0  /  HAVG 

C  =  0.5DOXBETAXEOEXSIN(XM(3))XSIN(XM(3))X(RE+HP) 

CC  =  C  X  C 

ZZ  =  BETA  X  XM(1)  X  XM(2) 

ALP  =  0.00 

MO  =  1 

NN  =  8 

CALL  MMBSIR(ZZ,ALP,NN,MO,BF,NZ) 

10  =  BF(1) 

11  =  BF(2) 

12  =  BF(3) 

13  =  BF(4) 
I<+  =  BF(5) 

15  =  BF(6) 

16  =  BF(7) 

17  =  BF(8) 
FA  =  FXA 
FB  =  F«B 

E  =  XM(2) 

C2W  =  C0S(2.XXM(4)) 

S2W  =  SIN(2.XXM(<+)) 

C4W  =  C0S(4.XXM(4)) 

S4W  =  SINC+.xXMCi)) 

DELMG=-DELX(SQRT(BOXXM(1)))XRHOOXEXP(-ZZ-CXC2W)X((1.+0.25DOX 
ICxC^x 

l(IO+EXIl)+FAx((l.+0.25DOXCC)x(li+Exl2))+0.5DOXCX((2.xi2-Ex(Il-3.X 
213  ))  +  FAx(Il  +  I3+2.xExi4))xC2W-0.5D0xcxFDX((Il-I3)  +  2.XEX( 12-14) )x 
3  S2W  +0.125D0XCXCX((2.XI4-CX(3.XI5-5.XI5))+FAX((I3+I5)+EX(3.XI6-I2 
<4   )))XC4W  +  0.125D0XCXCXFBX((I5-I3)  +  EX(I2-4.XI4+3.XI6))XS4H) 

XTl  =  DELX(SQRT(BOXXM(1)))XRHOO 

XT2  =  DELMG  /  XTl 

XT3  =  EXP(-ZZ)x(IO  +  EXIl) 

XT'i  =  XT2  /  XT3 

DELMG  =  1000.  X  DELMG  /  9.8D0 
K...   PRINT  X,  ZZ,  RHOO,  XTl,  XT2,  XT3,  XT'*,  DELMG  ,  DEL 

RETURN 

END 

X 
X 

X  '  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

3<  xxxxxxxxxxxxxxxxxxxxxxx^xxxxxx 

X'  XX  XX 

X  XX     SUBROUTINE  JAC60     XX 

X  XX  XX 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 
X 

SUBROUTINE  JAC60  (D, DS, A, AS, BLA,Z, F107, RHO, SH  ) 

BL  =  BLAX3.14159D0  /  180. 

XP  =  -16.021D0  -  0.001985D0XZ  +  6.363D0XEXP(  -0.0026D0XZ) 

RQ  —  10  XXXP 

CSPSI  ='SIN(D)XSIN(DS)  +  COS( D)xCOS( DS)XC0S( A-AS-BL) 

F  =  (O.SDOxd.  +  CSPSI))XX3. 

RHOl  =  1.  +  0.19D0X(EXP(0.0055D0XZ)  -  1.9D0)  XF 
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RHO  =  O.OIDO  X  F107  x  RO  X  RHOl  X  (10.XX12) 

DRODZ  =  2.3026D0XR0X(-.001985D0  -  .  01654'4D0xEXP(-.  0026D0XZ)  ) 

DRIDZ  =  .0010^5DOXFXEXP( .0055D0XZ) 

SH  =  -(ROKRH01)/(ROXDR1DZ+RH01XDRODZ) 

RETURN 

END 

X 

X  XJ(XXXX30(XXX3(X9(X3(XXXXXXX>(XXKXXX 

X  XX  *  XX 

X  XX  SUBROUTINE   THRST  XX 

X  XX  XX 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X  xxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

X 
X 

SUBROUTINE  THRSTC  XM.  DXM,  IRC,  K,  PA,  DWOA  ) 

DIMENSION  XM(6)  ,  DXM(6),  IRV(IO),  DA(IO) 

COMMON  /  OADJ  /  NOA,  IRV,  DA,  SPI,  WT,  PSENS 

DATA  BO  /  398603. 254D0/,  PI/3 .  l'il59D0/ 

DO  10   I  =  1,6 

DXM(I)  =0.0 
10  CONTINUE 

DWOA  =0.0 

IF  (  IRC  .EQ.  0)  RETURN 

THT  =  XM(6)  +  2.XXM(2)XSIN(XM(6))  +  (5 ./A . )XXM(2)XXMC2)XSIN(2 . XXM 
1      (6)) 

R  =  XMCnxCl.  -  XM(Z)XXM(2))/(1.  +  XM(2)xC0S(THT) ) 

V2  =  BO  X  ((2./R)  -  (l./XM(l))) 

V  =  SQRT(V2) 

DVI  =  DA(K)XB0/(2.XXM(1)XXM(1)XV) 

DXMC1)=  DA(K) 

DXM(2)  =  2.X((C0S(THT)+XM(2))/V)XDVI 

DXM(^)  =  2.XSIN(THT)XDVI/(XM(2)XV) 

DXM(6)  =  (2.XPI/PA)-(1./(XM(1)X(SQRT(1.-XM(2)XXM(2)))XV))X2.X 
1    SIN(THT)X((XM(1)X(1.-XM(2)XXM(2))/XM(2))  +  RXXM(2nxDVI 

PA  =  3.XPAXXM(1)X(V/B0)XDVI 

DWOA    =    (WT/CSPIX9.8D0))X1000.XDVI 
.    IF    (PSENS    .EQ.    0)    THEN 

PRINT    X,    THT,R,V,IRV(K),DA(K),DVI,DXM(2),DXM(^),DXM(6),PA, 
1  DWOA 

END    IF 

RETURN 

END 
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APPENDIX  B 
SAMPLE  INPUT  AND  OUTPUT 


05000,001,0,1,0,025 

2.2,  .0000500,  1.0,  230. ,100. ,09100. 

100.,  100.,  2.0  ,<+3222.7382 

1  109<i7  U  90070  7S   060    A  0        1.0  014791  1015        0        0  ++ 

2  10947  44619.98716775  150.0000  095.0000  200.0000  .0100000  095.0000  ++ 

3  13947  15.71356734   .000783912   1.00112  -3.81308  -0.23959E-00     0  ++ 

4  10947  000000000-0  0         0         0         0  ++ 

5  10947  01.06000000  -701.8652  E-00  ++ 

cccc  x.x  xxxxxx.xxxxxxxxx   xxxxxx.xxxxxxxxx 

CD  1.0  000001.000000000  000003.500000000 

A  0.0  000000.000010000  000000.001000000 

D  0.0  000000.100000000  000010.000000000 

SPI  0.0  000200.000000000  003000.000000000 

W  0.0  000005.000000000  000200.000000000 

WT  0.0  000100.000000000  100000.000000000 

F107  0.0  000050.000000000  000300.000000000 

FBAR  0.0  000050.000000000  000300.000000000 

AKP  0.0  000000. OCOOOOOOO  000007.000000000 

TVE  0.0  038422.000300000  039422.003000000 

UJD  0.0  043223.000000000  044223.000000000 

ETU  0.0  000000.000000000  000179.0000000  00 

KO  0.0  000000. OOOOOCOOO  000  36  0.000000000 

GO  0.0  000000.000000000  000360.000000000 

E5  0.0  000000.010000000  000000.100000000 

XI  0.0  000022.500000000  000157.500000000 

RND  0.0  000000.000010000  0000  00.000000000 

ESD  0.0  000003  .000000000  000  0  00.30000000  0 

AO  0.0  000001.025940000  000001.16347  00  00 

ADOT  0.0  000000.000000000  000000.000000000 

1000,  0.1 
2030,  0.3 
3000,  0.5 

X  X 

K   ONLY  THE  ABOVE  IS  ACTUAL  INPUT  FILE  DATA.   AN  INPUT  ELEMENT  X 

X   DESCRIPTION  IS  PRESENTED  NEXT  TO  AID  THE  USER  IN  PLACING  VALUES      x 

X  X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

INPUT  FILE  FOR  PROPELLANT  LONGEVITY  MODEL  (VARIABLE  ORGANIZATION) 

lETRl,  IETR2.  NOA,  PSENS,  lET,  PSIZE 
CD,  A,  D,  SPI,  W,  WT, 
F107,  FBAR,  AKP,  TVE, 

1  SKIP  LINE  ONE  NAVSPASUR  FIVE  CARD  DATA 

2  8X.  UJD,  IX,  ETU,  IX,  HO,  IX,  GO,  IX,  ES,  IX,  XI 

3  20X,  RND,  19X,  ESD 

4  SKIP  LINE  FOUR  NAVSPASUR  FIVE  CARD  DATA 

5  3X,  AO,  IX,  ADOT 
IRVCl),  DA(1)    1 

IRV(2),  DA(2)    !     NOTE:  WA.^NING  OMIT  ALL  ORBIT  ADJUS'  DATA  IF 
IRV(3),  DA(3)    I  NO  ORBIT  ADJUSTS  ARE  SCHEDULED  (NOA=0) 

I 

I 
I 

I 
IRV(X),  DA(K)    I 

SENSITIVITY  ANALYSIS  TABLE 

CHAR  SNSAR     SNSMIN       SNSMAX 
A4    F3.1      F16.8        F16.8 
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INITIAL  INPUT  CONDITIONS 


ATMOSPHERIC  INPUT  PARAMETERS 

D  (RELATIVE  ATMOSPHERIC  ROTATION)" =      1.0000 

F107  (DECIMETRIC  SOLAR  FLUX)  =   100.00 

FBAR  (90  DAY  AVG.  F1G7)  =   100.00 

AKP  (GEOMAGNETIC  INDEX)  =   2.000 

TVE  (TIME  OF  VERNAL  EQUINOX  PASSAGE)  =   <+3222 .  73S200000  MODIFIED  JULIAN  DAYS 


SATELLITE  PHYSICAL  PARAMETERS 

CD  (DRAG  COEFFICIENT)  =   2.20  UNITLESS 

A  (CROSSECTIONAL  AREA)  =  0.000050000  SQUARE  KILOMETERS 

SPI  (MOTOR  SPECIFIC  IMPULSE)  =    230.00  SECONDS 

W  (INITIAL  FUEL  MASS)  =      10.00  KILOGRAMS 

WT  (INITIAL  SATELLITE  MASS)  =      9100.00  KILOGRAMS 


NAVSPASUR  ORBITAL  ELEMENTS 

UJD  (EPOCH)  =  <i4619. 98716775  MEAN  JULIAN  DAYS 

ETU  (MEAN  ANOMALY)  =  150.0000  DEGREES 

HO  (MEAN  RIGHT  ASCENSION  OF  THE  ASCENDING  NODE)  =   95.0000  DEGREES 

GO  (MEAN  ARGUMENT  OF  PERIGEE)  =  200.0000  DEGREES 

XI  (MEAN  INCLINATION)  =   95.0000  DEGREES 

RND  (RATE  OF  CHANGE  OF  MEAN  MOTION)  =  0.000783912  REVOLUTTIONS  PER  DAY  SQUARED 

ES  (MEAN  ECCENTRICITY)  =    0.0100  UNITLESS 

ESD  (ECCENTRICITY  DECAY  RATE)  =      -0.239590000  1/SECONDS 

AO  (KAULA  SEMI-MAJOR  AXIS)  =      1.06000  MEAN  EARTH  RADII 

ADQT  (SEMI-MAJOR  AXIS  DECAY  RATE)  =    -701.865200000 


CONTROL  INPUT  PARAMETERS 

MAXIMUM  REVOLUTIONS  PER  ITERATION     5000 

NUMBER  OF  SCHEDULED  ORBIT  ADJUSTS  TO  SEMI-MAJOR  AXIS 


SENSITIVITY  ANALYSIS  FOR 

LOW  EARTH  ORBIT  SATELLITES 

PROPELLENT  LONGEVITY 


156 


TABLE  OF  INPUT  PARAMETERS  SCANNED  FROM  MINIMUM  TO 
MAXIMUM  WITH  THE  INDICATED  INCREMENT  PER  RUN. 
THE  GRANULARITY  IS  1/   25. 

ALL  OTHER  INPUT  PARAMETERS   ARE  RETURNED  TO  THEIR 
INITIAL  VALUES  AT  EACH  INCREMENT. 

INPUT  PARAMETER        MINIMUM  MAXIMUM  .  INCREMENT 

CD  1.000000000        3.500000000        O.IOOOOOOQO 


SENSITIVITY  SCAN  RESULTS  FOR  ABOVE  SPECIFICS 


REV  NUMBER 

TIME(DAYS) 

495 

31.58 

448 

28.57 

409 

26.08 

376 

23.97 

348 

22.18 

324 

20.65 

303 

19.30 

285 

18.15 

268 

17.07 

254 

16.17 

241 

15.34 

229 

14.57 

218 

13.87 

208 

13.23 

200 

12.72 

191 

12.14 

184 

11.70 

177 

11.25 

170 

10.30 

154 

10.42 

159 

10.10 

154 

9.78 

149 

9.46 

144 

9.14 

140 

8.88 

136 

8.63 

FRACTION  OF  SCAN 
.  1.000000000 
1.100000000 
1.200000000 
1.300000000 
>. 400000000 
1.500000000 
1.600000000 
1.700000000 
1.800000000 
1.900000000 
2.000000000 
2.100000000 
2.200000000 
2.300000000 
2.400000000 
2.500000000 
2.600000000 
2.700000000 
2.800000000 
2.900000000 
3.000000000 
3.100000000 
3.200000000 
3.300000000 
3.400000000 
3.500000000 
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05000,001,0,1,0,025 

2.2,  .0000500,  1.0,  250. ,001. ,09100. 

100.,  100.,  2.0  ,<+3222.7382 

1  10947  U  90070  78   060    A  0        1.0  014791  1015        0        0  ++ 

2  10947  44619.98716775  150.0000  095.0000  200.0000  .0100000  095.0000  ++ 

3  10947  15.71356734   .000733912   1.00112  -3.81308  -0.23959E-00     0  ++ 

4  10947  000000000-0  0         0         0         0  ++ 

5  10947  01.06000000  -701.8652  E-00  ++ 

cccc  x.x  xxxxxx.xxxxxxxxx        xxxxxx.xxxxxxxxx 

CD  0.0    000001.000000000  000003.500000000 

A  0.0    000000.000010000  000000.001000000 

D  0.0    000000.100000000  000010.000000000 

SPI  0.0    000200.000000000  003000.000000000 

W  0.0    000005.000000000  000200.000000000 

WT  0.0    000100.000000000  100000.000000000 

F107  0.0    000050.000000000  000300.000000000 

FBAR  0.0    000050.000000000  000300.000000000 

AKP  0.0    000000.000000000  00U007 . 000000000 

TVE  0.0    038422.000000000  039422.000000000 

UJD  0.0    043223.000000000  044223.000000000 

ETU  0.0    000000.000000000  000179.000000000 

HO  0.0    000000.000000000  000360.000000000 

GO  0.0    000000.000000000  000360.000000000 

ES  0.0    000000.010000000  000000.100000000 

XI  0.0    000022.500000000  000157.500000000 

RND  0.0    000000.000010000  000000.000000000 

ESO  0.0    000000.000000000  000000.000000000 

AO  1.0    000001.025940000  000001.168470000 

AOOT  0.0    000000.000000000  000000.000000000 

1000,     0.1 

2000,    0.3  ^ 

3000,    0.5 

X  ""  X 

X  ONLY   THE   ABOVE    IS    ACTUAL    INPUT    FILE    DATA.       AN    INPUT    ELEMENT  x 

X  DESCRIPTION    IS    PRESENTED   NEXT   TO    AID    THE   USER    IN    PLACING   VALUES  X 

X  X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

INPUT    FILE    FOR   PROPELLANT    LONGEVITY   MODEL    (VARIABLE   ORGANIZATION) 

lETRl,    IETR2,    NOA,    PSENS,    lET,    PSIZE 
CD,    A,    D,    SPI,    W,    WT, 
F107,    FBAR,    AKP,    TVE, 

1  SKIP    LINE   ONE   NAVSPASUR    FIVE   CARD    DATA 

2  3X,    UJD,    IX,    ETU,    IX,    HO,    IX,    GO,    IX,    ES,    IX,    XI 

3  20X,    RND,    19X,    ESD 

4  SKIP  LINE  FOUR  NAVSPASUR  FIVE  CARD  DATA 

5  8X,  AO,  IX,  ADOT 
IRV(l),  DA(1) 

IRV(2),  DA(2)    I     NOTE:  WARNING  OMIT  ALL  ORBIT  ADJUST  DATA  IF 
IRV<3),  DA(3)    I  NO  ORBIT  ADJUSTS  ARE  SCHEDULED  (NOA=0) 


IRV(K),  DA(K) 

SENSITIVITY  ANALYSIS  TABLE 

CHAR  SNSAR     SNSMIN       SNSMAX 
A4    F3.1      F16.8        F16.8 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 
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INITIAL  INPUT  CONDITIONS 


ATMOSPHERIC  INPUT  PARAMETERS 

D  (RELATIVE  ATMOSPHERIC  ROTATION)  =      1,0000 

F107  (DECIMETRIC  SOLAR  FLUX)  =   100.00 

FBAR  (90  DAY  AVG.  F107)  =   100.00 

AKP  (GEOMAGNETIC  INDEX)  =   2.000 

TVE  (TIME  OF  VERNAL  EQUINOX  PASSAGE)  =   43222 . 73S200000  MODIFIED  JULIAN  DAYS 


SATELLITE  PHYSICAL  PARAMETERS 

CD  (DRAG  COEFFICIENT)  =   2.20  UNITLESS 

A  (CROSSECTIONAL  AREA)  =  0.000050000  SQUARE  KILOMETERS 

SPI  (MOTOR  SPECIFIC  IMPULSE)  =    230.00  SECONDS 

W  (INITIAL  FUEL  MASS)  =       1.00  KILOGRAMS 

WT  (INITIAL  SATELLITE  MASS)  =      9100.00  KILOGRAMS 


NAVSPASUR  ORBITAL  ELEMENTS 

UJD  (EPOCH)  =  <44619. 98716775  MEAN  JULIAN  DAYS 

ETU  (MEAN  ANOMALY)  =  150.0000  DEGREES 

HO  (MEAN  RIGHT  ASCENSION  OF  THE  ASCENDING  NODE)  =   95.0000  DEGREES 

GO  (MEAN  ARGUMENT  OF  PERIGEE)  =  200.0000  DEGREES 

XI  (MEAN  INCLINATION)  =   95.0000  DEGREES 

RMD  (RATE  OF  CHANGE  OF  MEAN  MOTION)  =  0.000783912  REVOLUTTIONS  PER  DAY  SQUARED 

ES  (MEAN  ECCENTRICITY)  =    0.0100  UNITLESS 

ESD  (ECCENTRICITY  DECAY  RATE)  =      -0.239590000  1/SECONDS 

AO  (KAULA  SEMI-MAJOR  AXIS)  =      1.06000  MEAN  EARTH  RADII 

ADOT  (SEMI-MAJOR  AXIS  DECAY  RATE)  =    -701.365200000 


CONTROL  INPUT  PARAMETERS 

MAXIMUM  REVOLUTIONS  PER  ITERATION     5000 

NUMBER  OF  SCHEDULED  ORBIT  ADJUSTS  TO  SEMI-MAJOR  AXIS 


SENSITIVITY  ANALYSIS  FOR 

LOH  EARTH  ORBIT  SATELLITES 

PROPELLENT  LONGEVITY 
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TABLE  OF  INPUT  PARAMETERS  SCANNED  FROM  MINIMUM  TO 
MAXIMUM  HITH  THE  INDICATED  INCREMENT  PER  RUN. 
THE  GRANULARITY  IS  1/   25. 

ALL  OTHER  INPUT  PARAMETERS   ARE  RETURNED  TO  THEIR 
INITIAL  VALUES  AT  EACH  INCREMENT. 


INPUT  PARAMETER 
AO 


MINIMUM 
1.025940000 


MAXIMUM 
1.168470000 


INCREMENT 
0.005701200 


SENSITIVITY  SCAN  RESULTS  FOR  ABOVE  SPECIFICS 


REV  NUMBER 

TIME(DAYS) 

1 

0.00 

1 

0.00 

1 

0.00 

3 

0.12 

6 

0.31 

12 

0.70 

22 

1.34 

40 

2.51 

69 

4.42 

115 

,       7  .  47 

186 

12.21 

298 

19.76 

473 

31.65 

747 

■  50.42 

1149 

78.19 

1652 

113.33 

2033 

140.55 

2270 

158.15 

2514 

176.50 

2942 

208.13 

3349 

238.73 

3643 

262.00 

4026 

291.32 

4373 

318.80 

4555 

334.53 

4712 

348.62 

FRACTION  OF  SCAN 


025940000 
031641200 
037342400 
043043600 
048744800 
054446000 
060147200 
065848400 
1.071549600 
1.077250800 
1.082952000 
1.083653200 
1.094354400 


1.100055600 


105756300 
111453000 
117159200 
122360400 
1285616C0 
134262300 
1.139964000 
1.145665200 
1.151366400 
1.157067600 
1.162768800 
1.163470000 
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0500 
2.2, 

100. 
109 
109 
109 
109 
109 

CCCC 

CD 

A 

D 

SPI 

H 

WT 

F107 

F3AR 

AKP 

TVE 

UJD 

ETU 

HO 

GO 

ES 

XI 

RNO 

ESD 

AO 

ADQT 


0, 


001,0,1,0,025 

0000500,  1.0,  230. ,100. ,09100. 

100.,  2.0  ,^3222.7382 
U  90070  78   060    A  0        1.0  01^791  1015        0        0  ++ 
<+4619. 93716775  150.0000  095.0000  200.0000  .0100000  095.0000  ++ 
15.7135673A   .000783912   1.00112  -3.81308  -0.23959E-00     0  ++ 
000000000-0  0         0         0         0  ++ 

01.06000000  -701.8652  E-00  ++ 

X  xxxxxx.xxxxxxxxx   xxxxxx.xxxxxxxxx 

0  000001.000000000  000003.500000000 

0  000000.000010000  000000.001000000 

0  000000.100000000  000010.000000000 

0  000200.000000000  003000.000000000 

0  000005.000000000  000200.000000000 

0  000100.000000000  100000.000000000 

0  000050.000000000  000300.000000000 

0  000050.000000000  000300.000000000 

0  000000.000000000  000007.000000000 

0  038422. OOOOOOOCO  031422.000000000 

0  043223.000000000  044223.000000000 

0  000000.000000000  000179.000000000 

0  000000.000000000  000360.000000000 

0  000000.000000000  000360.000000000 

0  000000.010000000  000000.100000000 

0  000022.500000000  000157.500000000 

0  000000.000010000  000000.000000000 

0  000000.000000000  000000.000000000 

0  000001.025940000  000001.16847  0000 

0  000000.000000000  000000.000000000 


1000,  0.1 
2000,  0.3 
3000,  0.5 


X   ONLY  THE  ABOVE  IS  ACTUAL  INPUT  FILE  DATA.   AN  INPUT  ELEMENT  X 

X   DESCRIPTION  IS  PRESENTED  NEXT  TO  AID  THE  USER  IN  PLACING  VALUES       X 
X  X 

xxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxxx 

INPUT  FILE  FOR  PROPELLANT  LONGEVITY  MODEL  (VARIABLE  ORGANIZATION) 

lETRl,  IETR2,  NOA,  PSENS,  lET,  PSIZE 
CD,  A,  D,  SPI,  W,  WT, 
F107,  FEAR,  AKP,  TVE, 

1  SKIP  LINE  ONE  NAVSPASUR  FIVE  CARD  DATA 

2  3X,  UJD,  IX,  ETU,  IX,  HO,  IX,  GO,  IX,  ES,  IX,  XI 

3  20X,  RND,  19X,  ESD 

4  SKIP  LINE  FOUR  NAVSPASUR  FIVE  CARD  DATA 

5  SX,  AO,  IX,  ADOT 
IRV(l),  DA(1) 

NOTE:  WARNING  OMIT  ALL  ORBIT  ADJUST  DATA  IF 
NO  ORBIT  ADJUSTS  ARE  SCHEDULED  (NOA=0) 


IRV(2),  DA(2) 
IRV(3),  DA(3) 


IRV(K),  DA(K) 

SENSITIVITY  ANALYSIS  TABLE 


CHAR  SNSAR 
A4    F3 . 1 


SNSMIN 
F16.8 


SNSMAX 
F16.8 


XXXXXXXXXXXXXXXXXXXXXJ^XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX 
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INITIAL  INPUT  CONDITIONS 


ATMOSPHERIC  INPUT  PARAMETERS 

D  (RELATIVE  ATMOSPHERIC  ROTATION)  =      1.0000 

F107  (DECIMETRIC  SOLAR  FLUX)  =   100.00 

FBAR  (90  DAY  AVG .  F107)  =   100.00 

AKP  (GEOMAGNETIC  INDEX)  =   2.000 

TVE  (TIME  OF  VERNAL  EQUINOX  PASSAGE)  =   ^3222.738200000  MODIFIED  JULIAN  DAYS 


SATELLITE  PHYSICAL  PARAMETERS 

CD  (DRAG  COEFFICIENT)  =   2.50  UNITLESS 

A  (CROSSECTIONAL  AREA)  =  0.000050000  SQUARE  KILOMETERS 

SPI  (MOTOR  SPECIFIC  IMPULSE)  =    230.00  SECONDS 

H  (INITIAL  FUEL  MASS)  =      10.00  KILOGRAMS 

WT  (INITIAL  SATELLITE  MASS)  =     9100.00  KILOGRAMS 


NAVSPASUR  ORBITAL  ELEMENTS 

UJD  (EPOCH)  =  44619.98716775  MEAN  JULIAN  DAYS 

ETU  (MEAN  ANOMALY)  =  150.0000  DEGREES 

HO  (MEAN  RIGHT  ASCENSION  OF  THE  ASCENDING  NODE)  =   95.0000  DEGREES 

GO  (MEAN  ARGUMENT  OF  PERIGEE)  =  200.0000  DEGREES 

XI  (MEAN  INCLINATION)  =   95.0000  DEGREES 

RND  (RATE  OF  CHANGE  OF  MEAN  MOTION)  =  0.000783912  REVOLUTTIONS  PER  DAY  SQUARED 

ES  (MEAN  ECCENTRICITY)  =    0.0100  UNITLESS 

ESD  (ECCENTRICITY  DECAY  RATE)  =      -0.239590000  1/SECONDS 

AO  (KAULA  SEMI-MAJOR  AXIS)  =      1.06000  MEAN  EARTH  RADII 

ADOT  (SEMI-MAJOR  AXIS  DECAY  RATE)  =    -701.865200000 


CONTROL  INPUT  PARAMETERS 

MAXIMUM  REVOLUTIONS  PER  ITERATION     5000 

NUMBER  OF  SCHEDULED  ORBIT  ADJUSTS  TO  SEMI-MAJOR  AXIS 


SENSITIVITY  ANALYSIS  FOR 

LON  EARTH  ORBIT  SATELLITES 

PROPELLENT  LONGEVITY 
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TABLE   OF    INPUT    PARAMETERS    SCANNED    FROM   MINIMUM   TO 
MAXIMUM   WITH    THE    INDICATED    INCREMENT    PER   RUN. 
THE   GRANULARITY    IS    1/      25. 

ALL    OTHER    INPUT    PARAMETERS      ARE   RETURNED   TO   THEIR 
INITIAL    VALUES    AT    EACH    INCREMENT. 

INPUT    PARAMETER  MINIMUM  MAXIMUM  INCREMENT 

F107  50.000000000  300.000000000  10.000000000 


SENSITIVITY    SCAN    RESULTS    FOR   ABOVE   SPECIFICS 


REV   NUMBER 

TIME(DAYS) 

392 

24.99 

324 

20.65 

276 

17.58 

2«4l 

15.34 

213 

13.55 

191 

12.14 

174 

11.06 

159 

10.10 

146 

9.27 

136 

8.63 

127 

8.05 

119 

7.54 

112 

7.10 

105 

6.65 

100 

6.33 

95 

6.01 

90 

5.69 

86 

5.43 

82 

5.13 

79 

4.99 

75 

4.79 

73 

4.60 

70 

4.41 

67 

4.22 

65 

4.09 

63 

3.96 

FRACTION  OF  SCAN 
50.000000000 
60.000000000 
70.000000000 
80.000000000 
90.000000000 
100.000000000 
110.000000000 
120.000000000 
130.000000000 
140.000000000 
150.000000000 
160.000000000 
170.000000000 
180.000000000 
190.000000000 
200.000000000 
210.000000000 
220.000000000 
230.QOO0O0OOG 
240.000000000 
250.000000000 
260.000000000 
270.000000000 
280.000000000 
290.000000000 
300.000000000 
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INITIAL  INPUT  CONDITIONS 


ATMOSPHERIC  INPUT  PARAMETERS 

D  (RELATIVE  ATMOSPHERIC  ROTATION)  =     1.0000 

FI07  (DECIMETRIC  SOLAR  FLUX)  =   100.00 

F3AR  (90  DAY  AVG.  F107)  =   100.00 

AKP  (GEOMAGNETIC  INDEX)  =   2.000. 

TVE  (TIME  OF  VERNAL  EQUINOX  PASSAGE)  =   ^3222.738200000  MODIFIED  JULIAN  DAYS 


SATELLITE  PHYSICAL  PARAMETERS 

CD  (DRAG  COEFFICIENT)  =   2.20  UNITLESS 

A  (CROSSECTIONAL  AREA)  =  0.000050000  SQUARE  KILOMETERS 

SPI  (MOTOR  SPECIFIC  IMPULSE)  =    230.00  SECONDS 

W  (INITIAL  FUEL  MASS)  =     150.00  KILOGRAMS 

WT  (INITIAL  SATELLITE  MASS)  =     9100.00  KILOGRAMS 


NAVSPASUR  ORBITAL  ELEMENTS 

UJD  (EPOCH)  =  <+4619. 98716775  MEAN  JULIAN  DAYS 

ETU  (MEAN  ANOMALY)  =  150.0000  DEGREES 

HO  (MEAN  RIGHT  ASCENSION  OF  THE  ASCENDING  NODE)  =   95.0000  DEGREES 

GO  (MEAN  ARGUMENT  OF  PERIGEE)  =  200.0000  DEGREES 

XI  (MEAN  INCLINATION)  =   95.0000  DEGREES 

RND  (RATE  OF  CHANGE  OF  MEAN  MOTION)  =  0.000783912  REVOLUTTIONS  PER  DAY  SQUARED 

ES  (MEAN  ECCENTRICITY)  =    0.0100  UNITLESS 

ESD  (ECCENTRICITY  DECAY  RATE)  =     -0.239590000  1/SECONDS 

AO  (KAULA  SEMI-MAJOR  AXIS)  =      1.06000  MEAN  EARTH  RADII 

ADOT  (SEMI-MAJOR  AXIS  DECAY  RATE)  =    -701.865200000 


CONTROL  INPUT  PARAMETERS 

MAXIMUM  REVOLUTIONS  PER  ITERATION     5000 

NUM3ER  OF  SCHEDULED  ORBIT  ADJUSTS  TO  SEMI-MAJOR  AXIS 


PROPELLANT  LIFE  ESTIMATOR 
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CDA=  O.llOOOE-03  SQ  KM 

SPECIFIC  IMPULSE=  230.00  SEC 

INITIAL  WEIGHT  OF  PROPELLANT=  150.00  KG 


F107=   100.00 
FBAR=   100.00 
TIME  OF  VERNAL  EQaiNOX=   0 . 43222738200000E+05 


MAXIMUM  NUMBER  OF  ORBITAL  REVOLUTIONS=       5000 
DATA  PRINTED  EVERY  100  ORBITAL  REVOLUTIONS 

KOZAI  MEAN  ELEMENT  SET         "   • 

SEMIMAJOR  AXIS=   0  .  6756205'+107790E+0^  KM 
ECCENTRICITY=   0  .  96  5611329  0'+615E-02 
INCLINATION=   0  .  949999648<i9617E+02  DEG 
ARGUMENT  OF  PERIGEE=   0  .  193337'^2'i63643E+03  DEG 
ASCENDING  NODE=   0  .  95000149663'470E+02  DEG 
MEAN  ANOMALY=   0  .  156611526'*1309E+03  DEG 


REV  NUMBER 

TIMECMJD) 

TIME(DAYS) 

100 

44626.32 

6.33 

200 

44632.71 

12.72 

300 

44639.10 

19.11 

400 

44645.49 

25.50 

500 

44651.88 

31.90 

600 

44658.27 

38.29 

700 

44664.67 

44.68 

800 

44671.06 

51.07 

900 

44677.45 

57.46 

1000 

44683.84 

63.86 

1100 

44690.23 

70.25 

1200 

44696.63 

76.64 

1300 

44703.02 

33.03 

1400 

44709.41 

89.42 

1500 

44715.80 

95.82 

1600 

44722.19 

102.21 

1700 

44728.59 

108.60 

1300 

44734.93 

114.99 

1900 

44741.37 

121.33 

2000 

44747.76 

127.77 

2100 

44754.15 

134.17 

2200 

44760.55 

140.56 

2300 

44766.94 

146.95 

2400 

44773.33 

153.34 

2500 

44779.72 

159.73 

2600 

44786.11 

166.13 

.  2700 

44792.51 

172.52 

2800 

44798.90 

178.91 

2900 

44805.29 

185.30 

3000 

44811.68 

191.69 

3100 

44818.07 

198.09 

3200 

44324.47 

204.43 

3300 

44830.86 

210.87 

3400 

44337.25 

217.26 

3500 

44843.64 

223.65 

3600 

44850.03 

230.05 

3600 

230.05 

PROPELLANT  REMAINING(KG) 


145.28 

140.69 

136.21 

131.83 

127.54 

123.37 

119.33 

115.45 

111.72 

108.07 

104.42 

100.72 

96.94 

93.05 

88.96 

84.60 

79.92 

74.98 

69.93 

64.97 

60.27 

55.92 

51.91 

48.18 

44.59 

41.07 

37.58 

34.10 

30.53 

26.78 

22.74 

18.37 

13.66 

8. 68 

3.58 

-1.49 

0.000000000 
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